Application Serial No. 10/804,965 
Amendment and Response dtd. 1 1/2/07 
Response to Office Action of 8/6/07 

Remarks: 

The Office Action of August 6, 2007 has been carefully reviewed, and this response 
addresses the Examiner's concerns. 

I. STATUS OF THE CLAIMS 

Claims 1-24 are pending in the application. 

Claims 18-24 are rejected under 35 U.S.C. 101 as being directed to non-statutory subject 

matter. 

Claims 1-2, 4-9 and 1 1-24 are rejected under 35 U.S.C. 102(b) as being anticipated by 
Drisko et al. (US 5,933,523). 

Claims 3 and 10 are rejected under 35 U.S.C. 103(a) as being unpatentable over Drisko et 
al. and Yokoi (US 2001/0022854). 

II. OBJECTIONS TO THE SPECIFICATION 

Applicants enclose herewith a copy of J. S. Lim, Two Dimensional Signal and Image 
Processing, ISBN 0-13-935322-4, pp. 476-94 as requested by the Examiner. 

Applicants have also amended paragraph [0013] of the specification to remove the 
embedded hyperlinks. 

Finally, Applicants have amended the Abstract to be fewer than 150 words in length. 

III. THE 35 U.S.C. 101 REJECTION 

Claims 18-24 are rejected under 35 U.S.C. 101 as being directed to non-statutory subject 

matter. 

Applicants respectfully state that, for the reasons provided hereinbelow, claims 18-24 are 
patentable subject matter under 35 USC 101. 

In In re Beauregard (53 F. 3d 1583 (Fed. Cir. 1995)), the Commissioner stated that 
computer program products embodied in a computer usable medium are patentable subject 
matter under 35 USC 101. In re Beauregard followed the decision in In re Lowry (32 F. 3d 1 579 
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(Fed. Cir. 1994)) in which the Federal Circuit Court overturned a printed matter rejection for a 
claim for a memory including a data structure. 

In In re Lowry, The Board of Patent Appeals and Interferences held that the memory 
containing stored information recited an article of manufacture but did not give patentable weight 
to the data structure. The Court found that the data structures were "specific electrical magnetic 
structure elements in a memory and provided tangible benefits;" as such the data structures were 
physical entities. 

Following In re Lowry and In re Beauregard, a computer program product including a 
computer usable medium having computer readable code embodied therein is patentable subject 
matter and is covered by a statutory category. 

The use of the Internet for distribution and execution of software is a relatively new 
technology. See for example, Microsoft's decision to use the Internet for distribution of software 
as described in http://www.microsoft.com/presspass/press/1996/may96/esdpr.mspx or the use of 
CORBA or DCOM for distributed software execution involving a transmission from server to 
client which can happen via a carrier wave (the CORBA and DCOM definitions date to the mid- 
1990s) or the use of applets (also dating to the mid-1990s). 

The proposed interim guidelines stated that: 

from a technological standpoint, a signal encoded with functional descriptive 
material is similar to a computer-readable memory encoded with functional 
descriptive material, in that they both create a functional interrelationship with a 
computer. In other words, a computer is able to execute the encoded functions, 
regardless of whether the format is a disk or a signal. 

As such, a computer usable medium, having many possible embodiments, encoded with 
computer-readable code is a specific electrical/magnetic structure and provides a tangible result 
and, according to In re Lowry and In re Beauregard, is patentable subject matter and falls under a 
statutory category. 

Assuming, arguendo, that a computer usable medium as stated in In re Beauregard is 
only statutory for devices such as floppy disks, rigid magnetic disks, magnetic tape, optical tape, 
optical recording discs, punched tape or cards, or the like results in legal contradictions as 
presented herein below. 
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A. Taking arguendo the interpretation proposed by the Examiner, a computer usable memory 
has non-statutory equivalents under the Doctrine of Equivalents. (If a claimed invention is 
statutory, can the equivalents - that provide the same function, in the same way, to obtain the 
same result - be nonstatutory?) 

As the proposed guidelines state, a carrier wave having computer readable code embodied 
therein serves the same purpose, creating a functional interrelationship with a computer, as a 
computer readable memory and obtains the same result, the computer is able to execute the 
encoded functions. Also, a carrier wave having computer readable code embodied therein 
performs that same purpose in the same way as a computer readable memory having computer 
readable code embodied therein; both carry bits of information and both operate as 
communication channels (see for example, Wolf, J.K., Magnetic recording as a communications 
channel, 1994 IEEE International Symposium on Information Theory, 1994. Proceeding, 27 
June-1 July 1994 Page(s):5 or Moon, J, Signal-to-noise ratio definition for magnetic recording 
channels with transition noise, IEEE Transactions on Magnetics, Volume 36, Issue 5, Part 2, 
Sept 2000 Page(s):3881, 3883).' The carrier wave having computer readable code embodied 
therein obtains the same result, i.e. to enable a computer to provide a useful result by executing 
the code, as a computer readable memory having computer readable code embodied therein, 
An accused product or process may infringe a patent either literally or under the Doctrine of 
Equivalents. To literally infringe a patent claim, the product or process at issue must include 
each and every element of the claim. Builders Concrete, Inc v. Bremerton Concrete Products, 
757 F.2d 255, 257 (Fed. Cir. 1985). In the event of failing to meet the standards established for 
demonstrating a prima-facie case of literal infringement, the analysis then proceeds to another 
theory of infringement - the Doctrine of Equivalents. 

The Doctrine of Equivalents enables protect patent owners relief against infringers' 
making "unimportant and insubstantial changes" to a patented invention which "though adding 
nothing, would be enough ... [to evade] the reach of the law." Graver Tank & Mfg. Co. v.. Linde 
Air Products Co., 339 U.S. 605, 606, 70 S.Ct 854, 856, 94 L.Ed. 1097 (1950). The doctrine may 



1 Copies of the above cited pages of the Wolf and Moon articles are provided herewith in the Appendix and 

pertinent passages are highlighted. 
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result in a virtual expansion of the scope of a patentee's claims in certain circumstances where, 
although it is shown that the literal language of a claim's elements is not met by an accused 
product, it may be proven that the accused product, analyzed on an element by element basis, 
"performs substantially the same overall function or work, in substantially the same way, to 
obtain substantially the same overall result as the claimed invention." Pennwalt Co. v. 
Durand-Waylandlnc, 833 F.2d 931, 934 (Fed. Cir. 1987). 

Therefore, infringement under the Doctrine of Equivalents requires that each 
corresponding element of the accused apparatus, device, composition or process, that is not 
literally met by a limitation of a claim, performs substantially the same function in substantially 
the same way to yield the same, or substantially the same, result as each corresponding limitation 
of the claim. Graver Tank & Mfg. Co., Inc. v. Linde Air Products, Co., 339 U.S. 605; 
Perhin-Elmer Co. v. Computervision Corp. 732 F.2d 888, (Fed. Cir. 1984), cert, denied, 469 U.S. 
857 (1984); Pennwalt Co. v. Durand-Wayland, Inc., 833 F. 2d at 934 (Fed. Cir. 1987), cert, 
denied, 485 U.S. 961 (1988), and cert, denied, 485 U.S. 1009 (1988). 

The use of Doctrine of Equivalents can also be further burdened by file wrapper estoppel 
(Warner-Jenkinson Co. v. Hilton-Davis Chem. Co., 520 U.S. 17 (1997) and Festo Corp. v. 

Shoketsa Kinzoky Kogyo Kaboshiki Co. Ltd., 535 U.S. (2002); alleged disclaimers (Johnson 

& Johnson Assocs. v. U.R.E. Serv. Co., 285 F.3d 1046 (Fed. Cir. 2002) (en banc); and by prior 
art limits (Wilson Sporting Goods Co. v. David Geoffrey Assoc., 904 F.2d 677 (Fed. Cir. 1990) 
and if and when "means for" limitations per 35 U.S.C. § 1 12 (6 th par.) are deemed to be involved, 
by the interaction of that paragraph with the Doctrine of Equivalents as explained, e.g., in Al-Site 
Corp. v. VSI Intl Inc., 174 F.3d 1308 (Fed. Cir. 1999). 2 

As stated above and considering the law summarized here, a claim allowed before the 
mid-1990s reciting a computer usable memory having computer readable code embodied therein 
might protect against an infringing product utilizing a carrier wave (such as the Internet) to 
perform the same function since such a product would infringe under the Doctrine of Equivalents 



2 In Al-Site important distinctions are made between before-arising and after-arising technologies that can 

affect initial claim construction and can affect permissible range of application of the Doctrine of Equivalents. That 
distinction may also be significant in Doctrine of Equivalents consideration, even if no "means for..." limitation is 
involved. 



13 of 22 



Application Serial No. 10/804,965 
Amendment and Response dtd. 1 1/2/07 
Response to Office Action of 8/6/07 

unless subject to one or more of the other Doctrine of Equivalents limitations. Since the carrier 
wave is an after developed technology, limitations of the Doctrine of Equivalents place by 
Warner-Jenkinson and Festo and preceding cases might not apply to pre mid-1990s while those 
limitations of the Doctrine of Equivalents would apply to the Applicant's claims. 

Therefore, the application of the proposed guidelines to claims for computer program 
product comprising computer readable medium would render effective coverage of claims 
dependent on when the claim was written. This is an anomaly inconsistent with the clarifying 
purpose of the proposed guideline. Furthermore, the anti-carrier wave feature of the proposed 
guidelines presents a situation where, if the applicant does not explicitly list that a carrier wave is 
a computer usable medium, a claim for a computer usable medium, which is identical to the 
claims in In re Beauregard would protect against products using a carrier wave to embody 
computer readable code by means of the Doctrine of Equivalents while applicants who are more 
explicit in their definition of a computer usable medium would have their claims rejected, 
rendered patentable only by a narrowing amendment and therefore would be blocked by File 
Wrapper Estoppel from access to the Doctrine of Equivalents under Festo. 

The application of the Annex IV of the proposed guidelines to claims for computer 
program products comprising a computer usable medium presents a novel situation in patent law 
where a bona fide equivalent (under the Doctrine of Equivalents) is declared to be nonstatutory. 
Such a situation would be a case of first impression and is counter to applying the known law for 
determining compliance with 35 USC 101. 

B. The Cross Border Infringement Issue 

In NTP Inc. v. Research in Motion, Inc., 418 F.3d 1282, 75 USPQ 2d 1763 (Fed. Cir. 
2005) it was well explained that the steps of utilization of the Blackberry® system occurring 
partially in Canada did not infringe method claims of the patents in suit but that such usage did 
infringe apparatus (system) claims. 75 USPQ 2d at 1786-93. Nor did the defendant (nor 
customers of defendant induced by defendant) infringe by "exports" of text and voice message 
encoded signals to Canada under 35 U.S.C. § 271(f) [75 USPQ 2d at 1793-94] nor commit 
import of product-by-process infringement under 35 U.S.C. § 271(g) [75 USPQ 2d at 1794-95]. 
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Since under NTP v RIM methods claims do not provide protection against extra-territorial 

infringers, Applicant's protection would only be based on the system claim that is practiced by 

the user not the extra-territorial provider of the software. 

On the other hand in Eolas Tech Inc. v. Microsoft Corp., 399 F.3d 1375, 73 USPQ 2d 

1782 (Fed. Cir. 2005) and AT&T Corp. v. Microsoft Corp., 414 F.3d 1366, 75 USPQ 2d 1506 

(Fed. Cir. 2001) shipment of software on disks and indeed on a master disk and also as signals to 

Europe to be copied, with copies to be bundled with European made, sold, used computers did 

infringe a U.S. patent under the 35 U.S.C. § 271(f) export provisions. In AT&T, the Federal 

Circuit observed (75 USPQ 2d at 1509-10): 

"Were we to hold that Microsoft's supply by exportation of the master 
versions. . .avoids infringement we would be subverting the remedial nature of § 
271(f), permitting a technical avoidance of the statute by ignoring the advances in 
a field of technology - and its associated industry practices - that developed after 
the enactment of § 271(f). It would be unsound to construe a statutory provision 
that was originally enacted to encourage advances in technology by closing a 
loophole, in a manner that allows the very advances in technology thus 
encouraged to subvert that intent. Section 271(f), if it is to remain effective must 
therefore be interpreted in a manner that is appropriate to the nature of the 
technology at issue. . . ." 

The Federal Circuit's decision is consistent with a situation of streaming code abroad via 
the Internet rather than shipping a master disk. As stated by the Federal Circuit, "Additionally, 
we cannot accept Microsoft's suggestion that software sent by electronic transmission must be 
treated differently for purposes of § 271(f) liability from software shipped on disks, see Tr. of 
Dec. 12, 2003 Hearing, at 8:8-17 (J.A. 35 1), as it would amount to an exaltation of form over 
substance. "ATT v. Microsoft, 414 F. 3d 1366 (2005). 

A petition for Certiori was granted October 28, 2006 and Federal Circuit decision was 

reversed on other grounds. Microsoft v. ATT, 550 U. S. (2007), Decided April 30, 2007. As 

stated in the U.S. Supreme Court decision, "Until it is expressed as a computer-readable .copy,. 
e.g., on a CD-ROM, Windows software indeed any software detached from an activating 
medium remains uncombinable." Id. 

The patent in question in ,477/ v. Microsoft, RE32580, was re-issued in 1988 before the In 
re Lowry and in In re Beauregard decisions and therefore does not have a Beauregard type 
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claim. A Beauregard type claim would have avoided the controversy since, in claiming the 
software as being embodied in a computer readable medium (copy), the master disk (or the signal 
transmission, using the Federal Circuit statement) would have been the claimed invention and 
Microsoft would have been shipping an infringing product. 

Under the proposed Interim guidelines, transmission of software overseas to be bundled 
with a European made computer would not infringe while shipping a physical medium with the 
claimed software embodied therein would infringe. Therefore, the proposed Interim guidelines 
run counter to the statement by the Federal Circuit and permit "a technical avoidance of the 
statute by ignoring the advances in a field of technology." Also, in the words of the Federal 
Circuit, the proposed guidelines place to form over substance allowing infringers to use this form 
over substance guidelines to infringe without consequence. 

Should the Examiner apply guidelines not present in the MPEP to enable infringers to 
escape the statute by using advances in a field of technology that produce equivalents to the 
claimed invention? 

Applicants respectfully state that such a rejection would be counter to the Constitutional 
basis of the patent statute. 

For the reasons presented hereinabove, Applicants respectfully state that a utility rejection 
of claims 18-24 is not proper and should be withdrawn. 

IV. THE 35 U.S.C. 1 02(b) REJECTION 

Claims 1-2, 4-9 and 11-24 are rejected under 35 U.S.C. 102(b) as being anticipated by 
Driskoetal. (US 5,933, 523)(the '523 patent). 

Applicants respectfully assert that the claimed inventions are not anticipated by the '523 
patent since the '523 patent does not disclose culling the detected edges, detecting region corner 
points based on a predetermined relationship between each candidate corner point and 
characteristic edge points, and does not disclose obtaining a measure of cornerness (cornerness, 
also referred to as corner strength, is a term of art having a well known meaning for those skilled 
in the art) as detailed in the remarks given below. Regarding claim 1 1 , Applicants respectfully 
state that the limitations of claim 1 1 invoke 35 USC 112 paragraph 6; Applicants respectfully 
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state that the '523 patent does not disclose structures equivalent to the structures in the 
specification corresponding to the functions in the "means plus function" limitations of claim 1 1 . 

Regarding claims 1 and 18, the Examiner identifies the third limitation, "culling the 
detected edges in order to obtain a reduced edge group from the detected edges, the reduced edge 
group comprising a plurality of points," with Figures 4C and 4D and block 503 in Figure 5- 
Remove Known-Nozzle Points of the '523 patent. 

In the '523 patent, the operation at block 503 in Figure 5 is described in col. 9, lines 23- 
32 as follows: 

At 503, nozzle information is used to remove corner points corresponding to 
known nozzle points. Nozzle information may be in the form of a model, or a list 
of points in the physical or image space. Conventional techniques such as masking 
or point-by-point comparison may be used to identify and remove corner points 
corresponding to the nozzle at this stage. Similarly, information related to other 
objects known to show up in the image may be used to remove corner points 
believed to correspond to those objects. 

Applicants respectfully state that the operation block 503 in figure 5 up to '523 patent 
discloses removing corner points but does not disclose "culling the detected edges", which is a 
limitation of claim 1. 

The Examiner identifies "detecting region corner-points ...based on a predetermined 

relationship," another limitation of claim 1 , with block 502 of Figure 5 of the '523 patent. In the 

'523 patent, the extent coordinate system is described, in col. 7, lines 2739, as a linear 

transformation indicative of rotation, as shown below: 

At 501, the list of boundary points (x.sub.i, y.sub.i) are transformed from the 
image coordinate space into an extents coordinate space (x.sub.Ei, y.sub.Ei). The 
extents coordinate space is shown in FIG. 6C as a coordinate space which is fixed 
in relation to the coordinate space corresponding to the initial expected device 
orientation. Therefore, the axes (X.sub.E, Y.sub.E) of the extents coordinate space 
are offset by a fixed angle, for example, 45.degree., from the axes (X.sub.p, 
Y.sub.p) of the coordinate space corresponding to the image-space presentation 
angle (determined by mapping the physical space presentation angle via the 
calibration) representing the initial expected device orientation as it appears in the 
image. Because the relationship between the image coordinate space and the 
initial expected device orientation is described by the presentation angle and the 
relationship between the initial expected device orientation and the extents 
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coordinate space is fixed, a relationship between the image coordinate system and 
the extents coordinate system can be determined and described using common 
trigonometric formulas. 



In a linear transformation (in fact, in any one-to-one transformation whether linear or not) 
points in one space are mapped into correspondent points in the transformed space. Therefore, 
corner points in the image coordinate space are mapped into corresponding corner points in the 
extent coordinate space. 

As disclose in the '523 patent, in col. 8, lines 46-67 and in col. 9, lines 1-22 and in Figure 



6D: 



At 502, the boundary points transformed at 501 are evaluated to determine 
maximum (x.sub.max, y.sub.max) and minimum (x.sub.min, y.sub.min) values 
with respect to each of extents axes X.sub.E and Y.sub.E, which values define the 
extents. These extents are determined by minimum and maximum x and y 
positions of all the boundary points in the extents coordinate space. Both the 
extents and the extents points may be determined by comparing and updating the 
maximum and minimum values (x.sub.max, y.sub.max, x.sub.min, and y.sub.min) 
with each boundary point, x.sub.Ei, y.sub.Ei. The resulting extents points are 
shown, for example, as E1-E4 (in FIG. 6D). Extents points in the extents 
coordinate space correspond to corner points in the image coordinate space. 

Indices may be associated and updated with the extents points. Indices may be 
used to determine corner points within the image coordinate space corresponding 
to the extents points without requiring any further computation. For example, the 
following algorithm may be used to compute extents points and update 
corresponding indices: 

if Wto 

minX i=i 

else if x. iwx <x l .-, 

maxX_i=i 

if y mlt ,<y^ 

miitY i=i 

cLw tf y« B »<y ft v 
>•».„■..-» . 
nnixY i=i. 

In this algorithm, x.sub.min, x.sub.max, y.sub.min and y.sub.max correspond to 
extents in the extents coordinate space. The index maxX.sub.- i corresponds to 
the upper right corner of the chip, that given by minX.sub.-- i to the lower left 
corner, that given by maxY.sub.-- i to the upper left, and the that given by 
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minY.sub.-- i to the lower right. Using these relationships, corner points in the 
image coordinate space may be determined based on the extents points in the 
extents coordinate space without performing complex complete mapping or 
rotating operations. Alternatively, if the extents points are calculated using 
relationships (l)-(8) without using indexes, the extents points of the extents 
coordinate space may be translated into corner points of the image coordinate 
space using trigonometric relationships similar to relationships (l)-(8). 

In the above algorithm from the '523 patent, "extents" define the extent of the boundary 
of the device and corner points are mapped to corner points in the extent coordinate space and 
given indices so that corner points can be identified by the index without requiring the 
transformation. 

Applicants respectfully state that the 523 patent does not disclose "detecting region corner 
points from the plurality of candidate corner points based on a predetermined relationship 
between each of the candidate corner points and characteristic edge points," a limitation of claim 
1. 

Regarding claims 4 and 14, Applicants defined in the specification, in paragraph 13, 
cornerness as being the same as the term "corner strength" ("the measure of cornerness is also 
referred to as a measure of corner strength"). "Corner strength" is a term of art which has a well- 
defined meaning in the image processing/Computer vision field. (A number of references are 
provided in the Appendix showing widespread use of the term "corner strength" and in, some 
instances, the equivalent use of the term "cornerness.") Applicants respectfully state that the '523 
patent does not apply or disclose the use of the concept of corner strength as would be interpreted 
by one of ordinary skill in the art of image processing/computer vision. 

The first step in determining whether a claim is anticipated, or is obvious in view of prior 
art, is to interpret the claim. ("It is elementary in patent law that, in determining whether a patent 

is valid , the first step is to determine the meaning and scope of each claim in suit." 

Lemelson v. Gen. Mills, Inc., 968 F.2d 1202, 1206, 23 U.S.P.Q.2D (BNA) 1284, 1287 (Fed. Cir. 
1992).) When not defined by applicant in the specification, the words of a claim must be read as 
they would be interpreted by those of ordinary skill in the art. (MPEP 21 1 .01) (Rexnord Corp. v. 
Laitram Corp., 274 F.3d 1336, 1342, 60 USPQ2d 1851, 1854 (Fed. Cir. 2001) ("explaining the 
court's analytical process for determining the meaning of disputed claim terms")). Applicants 
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respectfully state that the term "cornerness" should be read as used in this specification and as 
would be understood by one skilled in the art, as evidenced by the references provided in the 
Appendix. 

The Examiner identified col. 6, lines 42-50 of the '523 patent as disclosing "obtaining a 
measure of cornerness," a limitation of claims 4 and 14. Col, 6, lines 42-50 of the '523 patent 
disclose: 

As described in more detail with respect to FIGS. 5 and 6A-6F, the illustrated 
embodiment specifically determines correspondence between the boundary 
features output at 302 and the corners of the generally rectangular device. For 
instance, as shown in FIG. 4D, the boundary features in FIG. 4C (output of the 
processing at 302) corresponding to the corners of the generally rectangular device 
are identified. As shown, only three of the boundary features (El, E2, E3) from 
the processing at 302 are identified as corresponding to the corners of the 
generally rectangular device since the boundary features most closely 
corresponding to the fourth corner actually correspond to the nozzle (depicted as 
402A in FIG. 4A). 

Applicants respectfully state that col. 6, lines 42-50 of the '523 patent do not disclose 
"obtaining a measure of cornerness" as would be understood by one of ordinary skill in the art. 

The Examiner identifies col. 6, line 5 of the '523 patent ("Once corresponding features 
are identified, they are fit against model input 414D to determine generally rectangular device 
position at 303B (same as 506 in FIG. 5), as shown in FIG. 4E.") as disclosing "selecting the 
plurality of candidate corner points from the plurality of point by applying a predetermined 
criterion," a limitation of claims 4 and 14. Applicants respectfully state that col. 6, line 5 of the 
'523 patent disclose determining generally rectangular device position which does not 
correspond to a claim limitation of claims 4 and 14. 

Regarding claims 5, 13 and 20, the comments provided above regarding the second 
limitation of claims 1 and 18 can be applied. Based on those comments, Applicants respectfully 
state that the 523 patent does not disclose the added limitations of claims 5, 13 and 20. 

Regarding claims 8, 16 and 23, based on the comments made above, Applicants 
respectfully state that the 523 patent does not disclose the added limitation of claims 8, 16 and 
23. 
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Regarding claim 1 1, the limitations of claim 1 1, 

means for detecting edges interior to a region of interest; 

means for culling the detected edges in order to obtain a reduced edge group from the 
detected edges, the reduced edge group comprising a plurality of points; 

means for selecting a plurality of candidate coiner points from the plurality of points; and, 

means for identifying region corner points from the plurality of candidate corner points 
based on a predetermined relationship between each candidate corner point and 
characteristic edge points of the region of interest; 

invoke 35 USC 1 12, paragraph 6. 

The MPEP in section 2182 details how to interpret (construe) a limitation that invokes 
35 U.S.C. 1 12, sixth paragraph. In order to construe a limitation invoking 35 U.S.C. 112, sixth 
paragraph, first, the function in the means plus function limitation has to be construed, and, 
second, the structures in the specification corresponding to the construed function have to be 
identified (MPEP 2182). 

Paragraph 13 and paragraph 15 of the Applicants' specification detail the structures 
corresponding to the recited functions in claim 1 1 . Applicants respectfully state, by the reasons 
given above, that the 523 patent does not disclose structure is equivalent to the corresponding 
structures in the Applicants' specification. 

"A claim is anticipated only if each and every element as set forth in the claim is found, 
either expressly or inherently described, in a single prior art reference." Verdegaal Bros. v. Union 
Oil Co. of California, 814 F.2d 628, 631, 2 USPQ2d 1051, 1053 (Fed. Cir. 1987); MPEP 2131. 

Claims 2, 4-9 are dependent on claim 1, claims 12-17 are dependent on claim 1 1 and 
claims 19-24 are dependent on claim 18. 

Applicants respectfully state that claims 1-2, 4-9 and 1 1-24 are not anticipated by the 
'523 patent. 
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V. THE 35 U.S.C. 103(a) REJECTION 

Claims 3 and 10 are rejected under 35 U.S.C. 103(a) as being unpatentable over Drisko 
et al. and Yokoi (US 2001/0022854)(the '854 publication). 

As stated above, the '523 patent (Drisko et al.) does not disclose all the limitations of 
claim 1. Claims 3 and 10 are dependent on claim 1. Applicants respectfully state that the Yokoi 
(the '854 publication) does not disclose the limitations of claim 1 that are not disclosed by the 
'523 patent. 

Since the '523 patent and Yokoi, either separately or in combination, do not teach or 
suggest all the limitations of claims 3 or 10, applicants respectfully state that a prima facie case 
of obviousness has not been established. 

VI. CONCLUSION 

In view of the above amendments and remarks, Applicants assert that claims 1-24 are 
now in condition for allowance and respectfully request the Examiner to pass the case to issue. 

The Director of Patents and Trademarks is authorized to charge any fee deficiencies, or to 
credit any overpayments, to Deposit Account No. 03-2410, Order No. 12078-205. 

In accordance with Section 714.01 of the M.P.E.P., the following information is 
presented in the event that a call may be deemed desirable by the Examiner: 



ORLANDO LOPEZ (617)345-3000 



Respectfully submitted, 

Lawrence E. Albertelli et al., Applicants 



Dated: November 2, 2007 




Orlando Lopez 
Reg. No. 46,880 
Attorney for Applicants 
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(a) 



(b) 



Figure 8.22 Example of salt-and-pepper noise reduction by median filtering, (a) Image in 
Figure 8.21(a) degraded by salt-and-pepper noise; (b) processed image by the same separable 
median filter used in Figure 8.21. 

resent additional details within the object may not be considered edges. An edge 
cannot be defined, then, outside of the context of an application. Nevertheless, 
edge detection algorithms that detect edges that are useful in a broad set of ap- 
plications have been developed. In this section, we discuss some of the more 
representative edge detection algorithms. 




Figure 8.23 Example of salt-and- 
pepper noise reduction by outrange 
pixel smoothing. Image in Figure 
8.22(a) processed by outrange pixel 
smoothing with threshold value of 50 
and window size of 3 x 3. 
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8.3.1 Gradient-Based Methods 



Consider an analog* function /(x) which represents a typical 1-D edge, as shown 
in Figure 8.24(a). In typical problems, it is reasonable to consider the value x Q in 
the figure an edge point. One way to determine x 0 is to compute the first derivative 
/' (*) or the second derivative f"{x). Figures 8.24(b) and (c) show /' (*) and /"(x) . 
From the figure, the value x 0 can be determined by looking for the local extremum 
(maximum or minimum) of /'(*) or by looking for a zero crossing of f"(x) where 
/"(*) changes its sign. In this section, we discuss methods that exploit the char- 
acteristics of /'(*). In the next section, we discuss methods that exploit the char- 
acteristics of f"(x). 

In addition to determining the possible edge point x 0 ,/'(x) can also be used 
in estimating the strength and direction of the edge. If |/'(jc)| is very large, f(x) 
is changing very rapidly and a rapid change in intensity is indicated. If /'(*) is 
positive, /(x) is increasing. Based on the above observations, one approach to 
detecting edges is to use the system shown in Figure 8.25. In the system, first 
|/'(*)| is computed from f(x). If [/'(*)) is greater than some threshold, it is a 
candidate to be an edge. If all values of x such that \f'(x)\ is greater than a certain 
threshold are detected to be edges, an edge will appear as a line rather than a 
point. To avoid this problem, we further require |/'(jc)| to have a local maximum 
at the edge points. It may be desirable to determine whether f(x) is increasing 
or decreasing at x = x 0 . The necessary information is contained in f'(x) at x = 
x 0 . The choice of the threshold depends on the application. As the threshold 
increases, only the values of x where f(x) changes rapidly will be registered as 
candidate edges. Since it is difficult to choose the threshold optimally, some trial 
and error is usually involved. It is also possible to choose the threshold adaptively. 
The system in Figure 8.25 is based on the particular type of edge shown in Figure 
8.24(a), but it is generally applicable to detecting various other types of edges. 

The generalization off (x) to a 2-D function f(x, y) is the gradient V/(* y) 
given by 

where \ x is the unit vector in the x-direction and t y is the unit vector in the y- 
direction. A generalization of the edge detection system in Figure 8.25 based on 
V/(x, y) is shown in Figure 8.26. The magnitude of Vf(x, y) is first computed 
and is then compared with a threshold to determine candidate edge points. If all 
values of (x, y) such that |V/(jc, y)| is greater than a certain threshold are detected 
to be edges, the edges will appear as strips rather than lines. The process of 
determining an edge line from a strip of candidate edge points is called edge 
thinning. In one simple edge thinning algorithm, the edge points are selected by 

*Sometimes, it is more convenient to develop results in the analog domain. In such 
instances, we will begin the development of results in the analog domain and then discretize 
the results at some later point in the development. 
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Figure 8.24 (a) /(*), (b) /'(*), and (c) f (x) for a typical 1-D edge. 



checking if |V/(jc, j)| is a local maximum in at least one direction. The property 
that |V/(x, _y)| achieves its local maximum in at least one direction is usually checked 
along a few specified directions. In most cases, it is sufficient to check for local 
maxima in only the horizontal and vertical directions. If |V/(x, y)\ is a local 
maximum along any one of the specified directions at a potential edge point, the 
potential edge point is considered to be an edge point. One difficulty with this 
simple edge thinning algorithm is that it creates a number of minor false edge lines 
in the vicinity of strong edge lines. One simple method to remove most of these 
minor false edge lines is to impose the following additional constraints: 



</<•] 


fix) 




i fix) 1 


> Threshold 


Yes 


Is 1 /'(x) ! a 


dx 















Figure 8.25 System for 1-D edge detection. 
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ix 0 , y 0 ) is not 
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Figure 8.26 System for 2-D edge detection. 



(a) If |V/(.r, y)| has a local maximum at (x 0 , y 0 ) in the horizontal direction but 
not in the vertical direction, (x 0 , y 0 ) is an edge point when 

k/(x,y)| Atffrt y)\ 

I dx I | dy J W1 tyP lcall y chosen around 2. 

(b) If |V/X*, y)\ has a local maximum at (x 0 , y 0 ) in the vertical direction but not 
in the horizontal direction, (x 0 , y 0 ) is an edge point when 

lal > H ^a'H with k ^P^^y ch °sen around 2. 



When |V/(jc, y)j has a local maximum at (x 0 , y 0 ) in the horizontal direction but not 
in the vertical direction, Condition (a) requires that the rate of intensity change 
along the horizontal direction is significantly larger than that along the vertical 
direction. Condition (b) is the same as Condition (a) with the roles of x and y 
reversed. Why these additional constraints remove most of the minor false edges 
in the vicinity of major edges is discussed in Problem 8.17. 

An edge detection system that is based on a function such as |V/(x, y)| is 
called a nondirectional edge detector, since such functions do not have a bias toward 
any particular direction. If an edge detection system is based on a function that 
has a bias toward one particular direction, it is called a directional edge detector. 
If we use \df(x, y)/dx\ instead of !V/(;t, y)| in the system in Figure 8.26, for example, 
the system will detect edges in the vertical direction, but will not respond to edges 
in the horizontal direction. 

For a 2-D sequence /(n ls n 2 ), the partial derivatives df(x, y)/dx and df(x, y)/dy 
can be replaced by some form of difference. For example, df(x, y)/dx may be 
replaced by 

~ [/(«!, «2) - - 1, n 2 )]/T, (8.7a) 
+ 1, «a) - n$/T, (8.7b) 
or [/(», + 1, n 2 ) - fin, - 1, ^{21). (8.7c) 

Since the computed derivatives are compared later with a threshold and the thresh- 
old can be adjusted, the scaling factors 1/7' and 1/2T can be omitted. Typically 
the expressions in (8.7) are averaged over several samples to improve the reliability 
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and continuity of the estimate of df(x, y)/dx. Examples of "improved" estimates 
of df(x, y)ldx are 



&U[/(« X + l,n 2 + I) -fin, - 1,« 2 + 1)] + [/(», + -fin, - l,n 2 )] 
+ [f( ni + 1, m - 1) - fin, - 1, « 2 - 1)] ( 8 - 8a ) 

or [/(«! + 1, » 2 + 1) - fin, ~l,n 2 + 1)] + 2[/(n, + 1,« 2 ) - /(% - l,«a)] 

+ [/(n, + 1, n 2 - 1) - /(»! - 1, « 2 - 1)1- t 8 " 8 ^ 

The unnecessary scaling factors have been omitted in (8.8). 

The differencing operation in (8.7) and (8.8) can be viewed as the convolution 
of fin,, n 2 ) with the impulse response of a filter h(n u n 2 ). Examples of impulse 
responses that can be used in developing directional edge detectors are shown in 
Figure 8.27. The filters kin,, n 2 ) in Figures 8.27(a) and (b) detect edges in the 
vertical and horizontal directions and can be viewed as approximating df(x, y)/dx 
and dfix, y)Idy respectively. The filters /i("i> "2) in Figures 8.27(c) and (d) detect 
edges along the two diagonal directions: The gradient V/(x, y) in (8.6) can also 
be expressed in terms of the first-order partial derivatives in a rotated coordinate 
system. When the rotation is 45 degrees, the directions of partial derivatives are 
along the two diagonal directions. 
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Figure 8.27 Impulse responses of fil- 
ters that can be used for directional 
edge detection, (a) Vertical edge detec- 
tion; (b) horizontal edge detection; 
(c) and (d) diagonal edge detection. 
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Nondirectional edge detectors can be developed by discrete approximation 
of |V/(at, y)\ in the system in Figure 8.26. From (8.6), 



From (8.9), nondirectional edge detectors can be developed by nonlinear combi- 
nation of the terms used in the development of directional edge detectors. An 
example of discrete approximation of (8.9) that can be used for nondirectional 
edge detectors is given by 

y)\ Vtffa, n 2 )f + (f y (n u n 2 )) 2 (8.10) 

where f x (n u n 2 ) = fin,, n 2 ) * h x (n u n 2 ) 

f y (n u n 2 ) = /(«!, n 2 ) * h y (n ly n 2 ) 

and h x (n lf n 2 ) and h y (n u n 2 ) are shown in Figure 8.28. The method developed 
by Sobel [Duda and Hart] is based on (8.10) with h x (n u n 2 ) and h y (n lt n 2 ) in Figure 
8.28. Another example is the method developed by [Roberts], which is based on 
(8.10) with h x (n u n^) and h y (n u n 2 ) shown in Figure 8.29. Depending on exactly 
how |V/(*, y)\ is approximated in the discrete domain, many other variations can 
be developed. 

Figure 8.30 shows the result of edge detection using a directional edge de- 
tector. Figure 8.30(a) shows an image of 512 x 512 pixels. Figure 8.30(b) and 
(c) show the results of a vertical edge detector and a horizontal edge detector, 
respectively, applied to the image in Figure 8.30(a). The vertical and horizontal 
edge detectors are based on h(n u n 2 ) in Figures 8.27(a) and (b). Figures 8.31(a) 
and (b) show the results of applying the Sobel edge detector and Roberts's edge 
detector to the image in Figure 8.30(a). Both belong to the class of nondirectional 
edge detectors, and the specific method of determining the threshold value and 
checking the local maximum property of an edge used is the same as that used in 
Figure 8.30. 
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Figure 8.28 Approximation of (a) 3f(x, y)tdx with /{n, , * , n 2 ) ; (b) df(x, y)l 
ty w ith /(fit, «z) * hyin-t, n 2 ). Sobel's edge detection meth od is based on comparison 
of V(/(n 1( nj) * kXn„ « 2 )) 2 + (/(«„ n 2 ) * h y {n„ n 2 )f with a threshold. 
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h y (n v n 2 ) 



Figure 8.29 Impulse responses of filters used in Roberts's edge detection method. 
The method is based on comparison of 

V(/(n„ n 2 ) * hAn„ n 2 )f + (/(„„ «J * h^n lt njf 
with a threshold. 

There are many variations of the edge detection methods discussed in this 
section . For example , we could use a different nonlinear combination of df(x, y)ldx 
and df(x, y)ldy instead of 

V(3/(*, y)/dxf + (a/(*, y)/dyf 

in the system in Figure 8.26. Many different methods can also be developed for 
edge thinning. 

The edge detection methods discussed in this section can be improved in 
various ways. Methods based on computing some form of gradient or differencing 
are typically sensitive to noise. A certain number of isolated edge points which 
appear randomly distributed throughout the edge maps in Figure 8.31 are most 
likely the result of some background noise or very fine image details. Some noise 
smoothing using the methods discussed in Section 8.2 or more sophisticated noise 
reduction methods that we will discuss in Chapter 9 may be desirable prior to 
applying an edge detection algorithm. Isolated random edge points may also be 
removed by some simple processing of the edge maps. Gradient-based edge de- 
tection methods also cause some discontinuities in the detected edge contours, as 
can be seen from the edge maps in Figure 8.31. Methods that impose continuity 
constraints in the detected edge contours can also be developed [Roberts]. 

8.3.2 Laplacian-Based Methods 



The objective of an edge detection algorithm is to locate the regions where the 
intensity is changing rapidly. In the case of a 1-D function f(x), searching for 
regions of rapidly changing intensity corresponds to searching for regions where 
f'(x) is large. For gradient-based methods, f'(x) is considered large when its 
magnitude is greater than a threshold. Another possible way is to assume 

that f'(x) is large whenever it reaches a local extremum, that is, whenever the 
second derivative f"(x) has a zero crossing. This is illustrated in Figure 8.24. 
Declaring zero-crossing points as edges results in a large number of points being 
declared to be edge points. Since there is no check on the magnitude of f'(x), 
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Figure 8.30 Edge maps obtained by directional edge detectors, (a) Image of 512 x 512 
pixels; (b) result of applying a vertical edge detector; (c) result of applying a horizontal edge 
detector. 
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Figure 8.31 Result of applying (a) Sobel edge detector and (b) Roberts's edge detector t< 
the image in Figure 8.30(a). 



any small ripple in/(jc) is enough to generate an edge point. Due to this sensitivity 
to noise, the application of a noise reduction system prior to edge detection is very 
desirable in processing images with background noise. 

A generalization of d 2 f(x)/dx 2 to a 2-D function f(x, y) for the purpose of 
edge detection (see Problem 8.19) is the Laplacian V 2 /(*, y) given by 

VY(*, y) = V(V/(, y) ) = ^> + mbJb. (8 . n) 

For a 2-D sequence f(n u n 2 ), the partial second derivatives d 2 f(x, y^dx 2 and 
d 2 f(x, y)Jdy 2 can be replaced by some form of second-order differences. Second- 
order differences can be represented by convolution of f(n u n 2 ) with the impulse 
response of a filter h(n u n 2 ). Examples of h(n t , n 2 ) that may be used are shown 
in Figure 8.32. To illustrate that f(n u n 2 ) * h(n u n 2 ) may be viewed as a discrete 
approximation of Wf(x, y), let us consider h(n u n 2 ) in Figure 8.32(a). Suppose 
we approximate df(x, y)tdx by 

-» fjn» » 2 ) = + 1, n 2 ) - /(«,, n 2 ). (8.12) 

We again omitted the scaling factor, since it does not affect zero-crossing points. 
Since the forward difference is used in (8.12), we can use the backward difference 
in approximating d 2 f(x, y^dx 2 : 

~* n 2 ) = f x (n u n z ) - fjfa - 1, n 2 ). (8.13) 
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Figure 8.32 Examples of n 2 ) that 
may be used in approximating V 2 f {x, y) 
with /(«„ « 2 ) * n 2 ). 



From (8.12) and (8.13), 
d 2 f(x, y) 

►/„(»!, « 2 ) = /(«, + 1, « 2 ) - 2/(n l5 « 2 ) + - l, „ 2 ). (8.14) 

From (8.11) and (8.14), and approximating d 2 f(x, y)/dy 2 in a similar manner, we 
obtain 

V^v)-*^/^,^) =/ i > lJ « 2 ) +/ w (« l3 n 2 ) 

= /(«! + 1, n 2 ) + /(rtj - 1, « 2 ) + f( ni , n z + 1) 

+ /(«i,«2-l)-4/(« 1 ,« 2 ). (8.15) 
The resulting V 2 f(n u n 2 ) is f(n u * h(n u n 2 ) with h(n u n 2 ) in Figure 8.32(a). 
Depending on how the second-order derivatives are approximated, it is possible 
to derive many other impulse responses h(n lt n 2 ), including those shown in Figures 
8.32(b) and (c). ' " 6 B 

Figure 8.33 shows an example where edges were detected by looking for zero- 
crossing points of V 2 /(r!, n 2 ). Figure 8.33(a) shows an image of 512 x 512 pixels. 
Figure 8.33(b) shows the zero-crossing points of V 2 /( ni , n 2 ), obtained from (8.15) 
and using the image in Figure 8.33(a) as f{n x , n 2 ). Since zero-crossing contours 
are boundaries between regions, they tend to be continuous lines. As a result, 
edge thinning necessary in gradient-based methods is not needed in Laplacian- 
based methods. In addition, algorithms that force continuity in edge contours are 
not as useful in Laplacian-based methods as in gradient-based methods. As is 
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clear from Figure 8.33(b), however, choosing ail zero-crossing points as edges tends 
to generate a large number of edge points. 

The Laplacian-based methods discussed above generate many "false" edge 
contours, which typically appear in regions where the local variance of the image 
is small. As a special case, consider a uniform background region so that/(n l5 n 2 ) 
is constant. Since V 2 /(n x , n 2 ) is zero and we detect edges from zero-crossing points 
of V 2 /("i, « 2 ), any small perturbation of f(n u n 2 ) is likely to cause false edge 
contours. One method to remove many of these false contours is to require that 
the local variance is sufficiently large at an edge point, as shown in Figure 8.34. 
The local variance ajin^ n 2 ) can be estimated by 

0 ^" 1 ' " z) = (2M + l) 2 M= i_M [ f(ku **> " m/{ku k ^ (816a) 

where mfa, n 2 ) = -^~f X f{h ' ^ {SMb) 

with M typically chosen around 2. Since oj(n 1; n 2 ) is compared with a threshold, 
the scaling factor 1/(2M + l) 2 in (8.16a) can be eliminated. In addition, the local 
variance aj needs to be computed only for (n u n 2 ) which are zero-crossing points 
of V 2 /(n 1; n 2 ). Figure 8.35 shows the result of applying the system in Figure 8.34 
to the image in Figure 8.33(a). Comparison of Figures 8.33(b) and 8.35 shows 
considerable reduction in the "false" edge contours. 

The system in Figure 8.34 can be interpreted as a gradient-based method. 




(a) (b) 



Figure 8.33 Edge map obtained by a Laplacian-based edge detector, (a) Image of 512 x 
512 pixels; (b) result of convolving the image in (a) with A(n„ n 2 ) in Figure 8.32(a) and then 
finding zero-crossing points. 
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i system that does not produce many false edge 



The local variance aj(n u n 2 ) is closely related to the gradient magnitude. Com- 
paring o}(n l5 n 2 ) with a threshold is similar to comparing the gradient magnitude 
with a threshold. Requiring that V 2 /(« a , n 2 ) crosses zero at an edge can be inter- 
preted as edge thinning. With this interpretation, we can implement the system 
in Figure 8.34 by computing oj(n lt n 2 ) first and then by detecting the zero-crossing 
points of V 2 /(n l7 n 2 ) only at those points where vjfa, ft 2 ) is above the chosen 
threshold. 



8.3.3 Edge Detection by Marr and Hildreth's Method 

In the previous two sections, we discussed edge detection algorithms that produce 
one edge map from an input image. Marr and Hildreth [Marr and Hildreth; Marr] 
observed that significant intensity changes occur at different scales (resolution) in 
an image. For example, blurry shadow regions and sharply focused fine-detail 




Figure 8.35 Edge map obtained by 
applying the system in Figure 8.34 to 
the image in Figure 8.33(a). 
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regions may be present in the same image. "Optimal" detection of significant 
intensity changes, therefore, generally requires the use of operators that respond 
to several different scales. Marr and Hildreth suggested that the original image 
be band-limited at several different cutoff frequencies and that an edge detection 
algorithm be applied to each of the images. The resulting edge maps have edges 
corresponding to different scales. 

Marr and Hildreth argue that edge maps of different scales contain important 
information about physically significant parameters. The visual world is made of 
elements such as contours, scratches, and shadows, which are highly localized at 
their own scale. This localization is also reflected in such physically important 
changes as reflectance change and illumination change . If the same edge is present 
in a set of edge maps of different scale, it represents the presence of an image 
intensity change due to a single physical phenomenon. ■ If an edge is present in 
only one edge map, one reason may be that two independent physical phenomena 
are operating to produce intensity changes in the same region of the image. 

To bandlimit an image at different cutoff frequencies, the impulse response 
h(x, y) and frequency response H(H X , fl y ) of the lowpass filter proposed [Marr 
and Hildreth; Canny] is Gaussian-shaped and is given by 

h(x,y) = e -(*»+j*WW) (g 17a ) 

H(£l x , ft,) = 2irVe-™ 3 < n * +sl *> /2 (8.17b) 
where cr determines the cutoff frequency with larger a corresponding to lower 
cutoff frequency. The choice of Gaussian shape is motivated by the fact that it is 
smooth and localized in both the spatial and frequency domains. A smooth h{x, y) 
is less likely to introduce any changes that are not present in the original shape. 
A more localized h{x, y) is less likely to shift the location of edges. 

From the smoothed images, edges can be detected by using the edge detection 
algorithms discussed in the previous two sections. Depending on which method 
is used, the lowpass filtering operation in (8.17) and the partial derivative operation 
used for edge detection may be combined. For example, noting that V 2 [-] and 
convolution * are linear, we obtain 

nf(x, y) * h(x, y)) = f{x, y) * [V*h(x, y)] 

(8.18) 



For the Gaussian function h(x, y) in (8.17), V 2 h(x, y) and its Fourier transform 
are given by 

V2/t (*> y) = (mr2)2 (* 2 + f - 2w 2 ) (8.19a) 
F[Wh(x, y)] = -Irf^e-^l+ma^ + a p (g J9b) 

Marr and Hildreth chose, for simplicity, to detect edges by looking for zero-crossing 
points of V 2 /(jc, y). Bandlimiting f(x, y) tends to reduce noise, thus reducing the 
noise sensitivity problem associated with detecting zero-crossing points. The func- 
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tions V 2 h(x, y) and -F[V 2 h(x, y)] in (8.19) are sketched in Figure 8.36. Clearly, 
computing /(x, y) * 7 2 h(x, y) is equivalent to bandpass filtering f(x, y) where o- 2 
in (8.19) is a parameter that controls the bandwidth of the bandpass filter. For a 
sequence f(n u n?), one approach is to simply replace * and y in (8.19) with n x and 

Figure 8.37 shows an example of the approach under discussion. Figures 
8.37(a), (b), and (c) show three images obtained by blurring the original image in 
Figure 8.33(a) with h(n u n 2 ) obtained by replacing x and y of h(x, y) in (8.17) with 
«! and n 2 with a 2 = 4, 16, and 36, respectively. Figures 8.37(d), (e), and (f ) show 
the images obtained by detecting zero crossings of/fa,, n 2 ) * V 2 fc(jc, y)\ x=nuy = n2 , 
with S7 2 h(x, y) given by (8.19a) for a 2 = 4, 16, and 36, respectively. Marr and 
Hildreth used the edge maps of different scales, such as those in Figures 8.37(d), 
(e), and (f ) for object representation in their image understanding work. 

8.3.4 Edge Detection Based on Signal Modeling 

The edge detection algorithms discussed above are general methods, in that they 
are developed independent of an application context. An alternative approach is 
to develop an edge detection algorithm specific to a particular application problem. 
If we know the shape of an edge, for example, this information can be incorporated 
in the development of an edge detection algorithm. To illustrate how an edge 
detection algorithm specific to an application problem may be developed, we con- 
sider the problem of detecting boundaries of coronary arteries from an angiogram 
[Abrams], 

The coronary arteries are the blood vessels which encircle the heart and supply 
blood to the heart muscle. Narrowing of the coronary arteries prevents adequate 
blood supply from reaching the heart, causing pain and damage to the heart muscle. 
Such damage is called coronary disease. To determine the severity of coronary 
disease , a coronary angiogram is used . An angiogram is an X ray picture of arteries 
taken after a contrast agent, typically iodine, has been injected into the vessels. 
The contrast agent is injected directly into the arteries through a catheter in order 
to achieve high concentrations. An example of a coronary angiogram is shown in 
Figure 8.38. Different observers making conventional visual evaluations of an 
angiogram will give widely varying evaluations of the severity of the disease. 

The most commonly used measure of an obstruction is percentage of stenosis, 
which is defined as the maximum percentage of arterial narrowing within a specified 
length of the vessel. One approach to estimating the percentage of stenosis begins 
with determining the vessel boundaries from an angiogram. We will be concerned 
with the problem of detecting the vessel boundaries. 

One reasonable model of an angiogram/^, n 2 ) is g iven by 

/(«!> «2> = (v(«i, « 2 ) + P(«i> * «2> + H"i> "2) (8-20) 

where v(n lt « 2 ) denotes the vessel, p(n u n 2 ) denotes the background, g(«x, rh) 
denotes blurring, and w(n u n 2 ) denotes the background noise. The vessel function 
v(n u n 2 ) is derived from a generalized cone model of a 3-D vessel which is con- 
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(a) (b) (el 




(d) (e) (f) 



Figure 8.37 Edge maps obtained from lowpass filtered image. Blurred image with (a) a 2 = 4; 
(b) o- 2 = 16; (c) a 2 = 36. Result of applying Laplacian-based algorithm to the blurred image; 
(d) a 2 = 4; (e) tr 2 = 16; (f) <r 2 = 36. 

tinuous and has elliptical cross sections. The elliptical shape is chosen because of 
the small number of parameters involved in its characterization and because of 
some empirical, evidence that it leads to a good estimate of percentage of stenosis. 
The 1-D cross section of vin^ n 2 ), which consists of one blood vessel, is totally 
specified by three parameters, two representing the blood vessel boundaries and 
one related to the x-ray attenuation coefficient of iodine. The continuity of the 
vessel is guaranteed by fitting a cubic spline function to the vessel boundaries. 
The background p(n u n 2 ) is modeled by a 2-D low-order polynomial. Low-order 
polynomials are very smooth functions, and their choice is motivated by the ob- 
servation that objects in the background, such as tissue and bone, are much bigger 
than the blood vessels. The blurring function g(nj, n 2 ) is modeled by a known 
2-D Gaussian function that takes into account the blurring introduced at various 
stages of the imaging process. The noise w^, n 2 ) is random background noise 
and is assumed to be white. The parameters in the model of f{n ly are the 
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Figure 8.38 Coronary angiogram. 



vessel parameters, the polynomial coefficients of p{n x , n 2 ), and the noise variance 
The vessels, tissues, bones, and the radiographic imaging process are much 
more complicated than suggested by the simple model presented above. Never- 
theless, the model has been empirically observed to lead to good estimates of the 
vessel boundaries and corresponding percentage of stenosis. The model param- 
eters may be estimated by a variety of different procedures. One possibility is the 
maximum likelihood (ML) parameter estimation method discussed in Section 6.1.5 
In the ML method, the unknown parameters denoted by e are estimated by max- 
imizing the probability density function PfMQ {f 0 { ni , n,)^) where /(«,, n 2 ) is 
the angiogram observation and 9 represents all the unknown parameters to be 
estimated. The ML method applied to vessel boundary detection is a nonlinear 
problem, but has been solved approximately [Pappas and Urn]. Figures 8.39 and 
8.40 illustrate the results of applying the ML parameter estimation method to the 
detection of the blood vessels using the 1-D version of the 2-D model in (8 20) 
In the 1-D version, f(n lt n 2 ) in (8.20) is considered a 1-D sequence with variable 
«i for each n 2 . Computations simplify considerably when the 1-D model is used 
Figure 8.39(a) shows the original angiogram of 80 x 80 pixels, and Figure 8 39(b) 
shows the detected vessel boundaries superimposed on the original image Figure 
8.40 is another example. Developing an edge detection algorithm specific to an 
application problem is considerably more complicated than applying the general 
edge detecdon algorithms discussed in previous sections. However, it has the 
potential of leading to much more accurate edge detection. 
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Magnetic Recording as a Communications Channel 

Jack Keil Wolf 

Center for Magnetic Recording Research, University or California, San Diego, La Jolla, CA 92093-0401, USA 

A communications channel is defined by specifying its in- 
put, its output a nd the probabilistic rules by which the al- 
lowable" inputs produce outputs. An oftejunsed, but overly 
simplistic, model ior a digital magnetic record ing channel is 
to descri be the input as a continuous waveform that takes on 
only" one of two val ues {siy +A arid -~A) and to describe the 
output as the sum of a linearly filtered versio n of the input 
and independent Gaussian noise. 



One aim of this talk is to explain how this simple model 
need be modified to take into account some of the important 
phenomena that occurs for high density recording. A second 
aim of this talk is to describe some of the coding schemes that 
are presently being used (01 are about to be used) in digital 
recording systems. 

My grandchildren sometimes axe rewarded with "special 
treats" when they are especially good. Your special treat for 
coming at this early hour on the last day of the conference (a 
schedule which was chosen by my good friend, the program 
committee chairman) will be to hear music played by the var- 
ious digital recording formats available today: i.e., CD, DCC, 
DAT and Minidisk. (This assumes that the necessary audio 
equipment will be available.) 
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Signal-to-Noise Ratio Definition for Magnetic Recording Channels with 
Transition Noise 

Jaekyun Moon, Senior Member, IEEE 



Abstract — This paper proposes a new signal-to-noise ratio 
(SNR) definition for magnetic recording channels with both 
additive and medium noise components. The proposed SNR is a 
generalized version of E h /N a , the information bit energy to noise 
spectral height ratio, widely used in average-power-constrained 
communication channels with additive white noise. The goal is 
to facilitate comparison of efficiencies of read channels that may 
operate at different symbol densities because of varying code rates. 



multiple errors, the detection SNR. of an optimal sequence 
detector is given by 



1 T J_ a 



[h(t)-h(t-T)] 2 dt 



Njzr 



(2) 



where T denotes the symbol period and E d is the dibit energy 
loss, medium noise, signal- 0 r the energy associated with a single nonreturn-to-zero (NRZ) 
written bit. This ratio can be interpreted as the ratio of the iso- 
lated dibit power to the in-band noise power (i.e., the noise 
I. Introduction power within the |/| < 1/2T band). It can be shown that E d 

decreases more rapidly than E s of an AWN, average power-con- 
strained channel, as the code rate decreases. This means that the 
penalty associated with a code rate loss is more severe with the 
magnetic recording channel than is an AWN channel with an 
average power constraint. For the Lorentzian channel, where 



We wish to define the signal-to-noise ratio (SNR) in such a 
way that we can compare efficiencies of different codes/detec- 
tors running at different symbol densities based on the SNR re- 
quired to meet a fixed bit error rate (BER) target. As such, the 
SNR definition should be free of the symbol density. However, 
the statistical properties of medium noise depend on the written 
pattern as well as on the symbol (or written bit) density, and 
defining an SNR to meet this goal is not a trivial matter. 

Let us first r eview SNR definitions for addi tive noise chan- 
nels. In communication channels with additive white noise 
(AWN) and an average transmitter power c onstraint, the E b /N„ 
ratio, the information bit energy to noise spectral height ratio, 
serves this purpose. It is well known that the error probability 
of such systems is mainly determined by the ratio 



~ N 0 



_ rE b _ 
= N 0 ~ 



P s 
N 0 /T' 



(1) 



where E 3 is the symbol energy, r is the code rate, and P s denotes 
the average transmitter power. This ratio is sometimes referred 
to as the detection SNR. The penalty associated with the code 
rate loss is apparent in (I). 

For magnetic recording channels with AWN, E t /N„, the iso- 
lated transi tion pulse energy to noise spectra l height ratio, can 
serve t his function as discussed in [1], A ma gnetic recording 
channel c an be viewed as the "^(t)-constrain ed ' channel, 
whe re h(t) is the channel's respon se to a single transition; there 
is no "average or peak power constraint here, but"for a given 
head/medium interface h(t) is fixed, and coding and signal 
processing engineers must live with it. At low-to-medium 
recording densities, where single errors are more likely than 
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h(t) = 



-I 



Vg 

1 -f- (2i/PW50) 2 



4Et 1 
tcPWSO ' 1 + (2t/PW50) 2 



(4) 



where D e denotes the symbol density defined as PW50/T and 
D u is the user density given by r ■ D s . It can be seen that /? 1Ilag is 
roughly proportional to r 2 (assuming D\ is considerably larger 
than r 2 ), as noted in [ 1 ], 1 whereas 0 com is proportional to r, con- 
firming a more severe code rate penalty in magnetic recording 
than in power constrained communication channels. 

II. Magnetic Recording Channels with Transition 
Noise 

The Et/N Q definition, however, cannot be used when transi- 
tion noise exists. We propose the following SNR definition to 
handle mixed noise channels: 
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where M 0 is twice the average energy 2 of the medium noise 
associated with each fransition, N n = No + Mn, and a denotes 

'Ryan used the notation to denote the energy in the isolated transition 
response, but we use the subscript t to emphasize it is the energy in the isolated 
transition response rather than isolated dibit response. Also, note thai h(i) here 
is the response to a full magnetization step, whereas the transition response in 
[I] is the half magnetization step, leading to the relationship: E t = 4Ei. 

2 It is safe to assume that the medium noise voltage associated with each tran- 
sition is square integrable. 
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the medium noise power expressed as a percentage of the total 
in-band noise power, i.e., 



M 0 
' No + Mo 



< ion. 



(6) 



We can further write No = [(100 - a)/100]N a and M 0 = 
(a/100)N a with 0 < a < 100. With a = 0, the SNR reduces to 
Et/No. Letting n m (f) denote the medium noise voltage wave- 
form associated with a given transition, the average medium 
noise energy per transition can be defined as 



Mo/2 = f" nl 



(t)dt, 



(7) 



where the statistical average is taken over transitions written at 
different positions in the disk. If the medium noise is statistically 
independent among transitions (as in the case of disks exhibiting 
a linear increase of noise power versus linear density), the total 
medium noise power originating from the all-transitions (or IT) 
pattern is easily shown to be nj*,(f) dt/T. This means that 



Ma 
2T = 



(total medium noise power in IT pattern) (8) 



and Mo can be viewed as the equivalent (single-sided) band-lim- 
ited white noise spectrum of the medium noise for this pattern 
(i.e., Mo is the height of a box-car power spectral density (PSD) 
of width 1/2T that integrates to the same value as the IT pat- 
tern's medium noise PSD). 

The SNR definition of (5) can be rewritten as 



SNR = 



E t _ 1 E t /T 
N n + M 0 ~ 2 " No Mo ' 
2T 2T 



giving rise to an interpretation 



SNR = 



Et/T 



2 (total inband noise power with IT pattern) 



(10) 



which can be viewed as the ratio of the isolated transition pulse 
power to the overall in-band noise power. 

For the first-order position jitter and width variation (1PW) 
noise model [2], [3], Mq can be written as 



Mo = 2aU t 4 



(U) 



where <r| and are the variances of position jitter and width 
variation, respectively, and 



/, = [dhfdtf dt, 

i m = j°° [dh/dw] 2 dt. 



(12) 
(13) 



With the above definitions, 5 specifying E t /N a with a partic- 
ular value of a determines both N a and M 0 [assuming h(t) is 
given as well]. Once M 0 is determined, the noise parameters 
necessary for BER simulation or analysis are specified. For ex- 
ample, with the 1 PW model, both a? and are specified once 

3 Here, we are implying h(t) should really be written as h{t, w), where w is 



Mo is known, provided the ratio of position jitter noise power 
to the width variation noise power is also known. Once N„, 
a-j, and, af„ as well as the symbol period are specified, we can 
proceed with performance analysis or BER simulations to com- 
pare codes/detectors operating at different symbol or user den- 
sities. With the microtrack noise model [4], fixing M 0 will de- 
termine the number of microtracks and the magnetization width 
parameter, which are required to simulate the recording channel 
output. Thus, assuming h(t) and the noise model are predeter- 
mined, specifying E t /N a with a particular value of a com- 
pletely characterizes the channel for analysis and simulation 
purposes. 

Note that the key assumption here is that M 0 is independent 
of the symbol period T so that the noise composition param- 
eter a is also independent of T. This means that the E t /N a 
definition is free of symbol density, as desired. This assumption 
will hold for any recording systems in which the total integrated 
medium noise power increases linearly with 1/T(i.e., where the 
intrinsic medium noise power per transition remains the same 
as the symbol pattern or density changes). This condition is met 
with well-designed media that do not exhibit supralinear noise 
increase with density (see, for example, [5]). 

For the Lorentzian channel, it can be shown that 4 



--Vo- 



(14) 



EJN„ 



(9) 
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'(15) 

As a point of reference, to meet a 10 ° target BER for the 1 PW 
Lorentzian channel at user density 2.5, EPR4ML (extended 
partial response class 4 maximum likelihood) with a rate 16/17 
runlength-limited (RLL) code requires an E t /N 0 slightly 
above 24 dB and an E t /Nso approximately equal to 23 dB with 
medium noise comprising only position (independent) jitter [6]. 

III. Relationships with Other SNR Definitions 
Let us compare the proposed E t /N a ratio with the SNR def- 
inition recording engineers are more familiar with, namely, the 
peak isolated pulse amplitude to the total in-band noise ratio for 
a pattern written at some particular density, say, D s . This SNR 
can be expressed as 

f^D.«) = jn £-^, (16) 

where V a is the peak-zero amplitude of h(t), NP(D,,-y) is 
the total l/2T-band noise power measured with a data pat- 
tern whose symbol density is D 3 , and average transition den- 
sity (total number of transitions divided by the total number of 
written bits) is 7. We can further write 



SNR P (£> B , 7 ) = 



y2 



E t 



7M0 
2T 



(17) 
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where the last equality holds for the Lorentzian pulse. Com- 
paring (17) and (5), the two SNRs can be related via 

Another SNR definition of interest is the mean-square (ms) SNR 
ratio defined as 



SNR 11)S (£> S ,7) = 




where A is some large interval (A ^> T), and .9(4) and n(t) 
are overall signal and noise (band-limited to 1/2T) voltages, re- 
spectively, corresponding to a written pattern with symbol den- 
sity D s and average transition density 7. For Lorentzian, we 
have 

™^D,n)- {D2s + l £ Mo + Noy (20) 

This SNR can be relate to E t /N a via 

SNR»(^,7 = 1) = ^ I ~- (2D 



IV. Conclusions 
A new SNR d efinition has been proposed tha t can be used to 
comp are read channels with different code rates. It can handle 
channels with medium noise. The proposed SNR definition, de- 
noted by Et/N a with a representing the percentage medium 
noise power, is applicable to any medium noise model as long as 
the model yields a linear relationship between the total medium 
noise power and the number of transitions in a given length of 
a track. The relationships between E t /N a and other often-used 
SNR ratios, namely, peak-zero SNR and ms SNR, have been es- 
tablished. 
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Abstract. There are many applications such as stereo matching, mo- 
tion tracking and image registration that require so called 'corners' to 
be detected across image sequences in a reliable manner. The Harris cor- 
ner detector is widely used for this purpose. However, the response from 
the Harris operator, and other corner operators, varies considerably with 
image contrast. This makes the setting of thresholds that are appropri- 
ate for extended image sequences difficult, if not impossible. This paper 
describes a new corner and edge detector developed from the phase con- 
gruency model of feature detection. The new operator uses the principal 
moments of the phase congruency information to determine corner and 
edge information. The resulting corner and edge operator is highly local- 
ized and has responses that are invariant to image contrast. This results 
in reliable feature detection under varying illumination conditions with 
fixed thresholds. An additional feature of the operator is that the corner 
map is a strict subset of the edge map. This facilitates the cooperative 
use of corner and edge information. 

1 Introduction 

With, the impressive reconstruction results that have been achieved by those 
working in projective geometry (see for example Hartley and Zisserman [1]) there 
has been a renewed interest in the detection of so called 'corners', or 'interest 
points'. The success of these reconstructions depend very much on the reliable 
and accurate detection of these points across image sequences. 

The definition of a corner is typically taken to be a location in the image 
where the local autocorrelation function has a distinct peak. A variety of op- 
erators have been devised to detect corners. These include those developed by 
Moravec [2], Harris and Stephens [3], Beaudet [4], Kitchen and Rosenfeld [5], 
and Cooper et al. [6]. Corner detectors based on the local energy model of fea- 
ture perception have been developed by Rosenthaler et al. [7], and Robbins and 
Owens [8]. More recently the SUSAN operator has been proposed by Smith and 
Brady [9] . Of these the Harris operator probably remains the most widely used. 

A common problem with all these operators, except the SUSAN operator, 
is that the corner response varies considerably with image contrast. This makes 
the setting of thresholds difficult. Typically we arc interested in tracking fea- 
tures over increasingly extended image sequences. The longer the sequence the 



greater the variations in illumination conditions one can expect, and the set- 
ting of appropriate thresholds becomes increasingly difficult, if not impossible. 
Inevitably thresholds have to be set at levels lower than ideal because detecting 
too many features is a lesser evil than detecting not enough. Stereo and motion 
reconstruction algorithms are then faced with the problem of dealing with very 
large clouds of noisy corner points, often greatly compromising their operation. 
Matching operations, which usually have rapidly increasing running time as a 
function of input size, suffer greatly in these conditions. Typically, considerable 
effort has to be devoted to the cleaning up of the output of the corner detector 
and to the elimination of outliers. The success of this cleaning up and outlier 
elimination process is usually crucial to the success of the reconstruction algo- 
rithm. Indeed, a very significant proportion of Hartley and Zisserman's book [1] 
is devoted to robust estimation techniques. 

Another difficulty many of these operators have is that the Gaussian smooth- 
ing that is employed to reduce the influence of noise can corrupt the location 
of corners, sometimes considerably. The SUSAN operator deserves some special 
comment here because it does not suffer from these problems outlined above. 
It identifies features by determining what fraction of a circular mask has values 
the same, or similar, to the value at the centre point. Thresholds are therefore 
defined in terms of the size of the mask and no image smoothing is required. 
However, the SUSAN operator assumes that edges and corners are formed by 
the junctions of regions having constant, or near constant, intensity, and this 
limits the junction types that can be modeled. 

To address the many problems outlined above this paper describes a new 
corner and edge detector developed from the phase congruency model of feature 
detection. The new operator uses the principal moments of the phase congruency 
information to determine corner and edge information. Phase congruency is a 
dimensionless quantity and provides information that is invariant to image con- 
trast. This allows the magnitudes of the principal moments of phase congruency 
to be used directly to determine the edge and corner strength. The minimum 
and maximum moments provide feature information in their own right; one does 
not have to look at their ratios. If the maximum moment of phase congruency 
at a point is large then that point should be marked as an edge. If the mini- 
mum moment of phase congruency is also large then that point should also be 
marked as a 'corner'. The hypothesis being that a large minimum moment of 
phase congruency indicates there is significant phase congruency in more than 
one orientation, making it a corner. 

The resulting corner and edge operator is highly localized and the invariance 
of the response to image contrast results in reliable feature detection under 
varying illumination conditions with fixed thresholds. An additional feature of 
the operator is that the corner map is a strict subset of the edge map. This 
facilitates the cooperative use of corner and edge information. 

This paper is organized as follows: first the phase congruency model of feature 
perception is reviewed. We then examine how the phase congruency responses 
over several orientations can be analyzed in terms of moments to provide both 



edge and corner information. Finally the performance is assessed relative to the 
commonly used Harris operator. 

2 The Phase Congruency Model of Feature Detection 

Rather than assume a feature is a point of maximal intensity gradient, the Local 
Energy Model postulates that features are perceived at points in an image where 
the Fourier components are maximally in phase as shown in Figure 1 [10]. Notice 




how the Fourier components are all in phase at the point of the step in the square 
wave. Congruency of phase at any angle produces a clearly perceived feature [11]. 
The angle at which the congruency occurs dictates the feature type, for example, 
step or delta. 

The Local Energy Model model was developed by Morrone et al. [10] and 
Morrone and Owens [12]. Other work on this model of feature perception can be 
found in Morrone and Burr [13], Owens et al. [14], Venkatesh and Owens [15], 
and Kovesi [16-20, 11]. The work of Morrone and Burr [13] has shown that this 
model successfully explains a number of psychophysical effects in human feature 
perception. 

The measurement of phase congruency at a point in a signal can be seen 
geometrically in Figure 2. The local, complex valued, Fourier components at a 
location x in the signal will each have an amplitude A n (x) and a phase angle 
<p n {x). Figure 2 plots these local Fourier components as complex vectors adding 
head to tail. The magnitude of the vector from the origin to the end point is the 
Local Energy, \E(x)\. 

The measure of phase congruency developed by Morrone et al. [10] is 

Under this definition phase congruency is the ratio of \E(x)\ to the overall path 
length taken by the local Fourier components in reaching the end point. If all 
the Fourier components are in phase all the complex vectors would be 
and the ratio of \E(x)\/ £ n A n {x) would be 1. If there is no coherence o 



Pig. 2. Polar diagram showing the Fourier components at^a location in the signal plotted 
head to tail. The weighted mean phase angle is given by <f>{x). The noise circle represents 
the level of E(x) one can expect just from the noise in the signal. 



the ratio falls to a minimum of 0. Phase congruency provides a measure that is 
independent of the overall magnitude of the signal making it invariant to vari- 
ations in image illumination and/or contrast. Fixed threshold values of feature 
significance can then be used over wide classes of images. 

It can be shown that this measure of phase congruency is a function of the 
cosine of the deviation of each phase component from the mean 



This measure of phase congruency does not provide good localization and it is 
also sensitive to noise. Kovesi [18, 19] developed a modified measure consisting of 
the cosine minus the magnitude of the sine of the phase deviation; this produces a 
more localized response. This new measure also incorporates noise compensation: 

PC2{X) E„^) + e • 

(3) 

The term W(x) is a factor that weights for frequency spread (congruency over 
many frequencies is more significant than congruency over a few frequencies). A 
small constant, e is incorporated to avoid division by zero. Only energy values 
that exceed T, the estimated noise influence, are counted in the result. The 
symbols [ J denote that the enclosed quantity is equal to itself when its value is 
positive, and zero otherwise. In practice local frequency information is obtained 
via banks of Gabor wavelets tuned to different spatial frequencies, rather than via 
the Fourier transform. The appropriate noise threshold, T is readily determined 
from the statistics of the filter responses to the image. For details of this phase 
congruency measure and its implementation see Kovesi [19-21]. 



3 Combining Phase Congruency Information over many 
Orientations 

A weakness of previous implementations of phase congruency has been the way 
in which information over many orientations is used and combined. The defini- 
tion of phase congruency outlined above only applies to ID signals. To obtain 
an overall measure of phase congruency in 2D local energy is first calculated in 
several orientations, typically six, using data from oriented 2D Gabor wavelets. 
Equation 3 is modified so that the numerator is the weighted and noise com- 
pensated local energy summed over all orientations, and the denominator is the 
total sum of filter response amplitudes over all orientations and scales. While 
this approach produces a phase congruency measure that results in a very good 
edge map it ignores information about the way phase congruency varies with 
orientation at each point in the image. 

To include information about the way phase congruency varies with orienta- 
tion we can proceed as follows: calculate phase congruency independently in each 
orientation using equation 3, compute moments of phase congruency and look at 
the variation of the moments with orientation. The principal axis, corresponding 
to the axis about which the moment is minimized, provides an indication of the 
orientation of the feature. The magnitude of the maximum moment, correspond- 
ing to the moment about an axis perpendicular to the principal axis, gives an 
indication of the significance of the feature. If the minimum moment is also large 
we have an indication that the feature point has a strong 2D component to it, 
and should therefore be additionally classified as a 'corner'. 

Following the classical moment analysis equations [22] we compute the fol- 



lowing at each point in the image: 

a=£(PC(0)cos(0)) 2 (4) 

6 = 2 £(PC?(0) cos(e)).(PC(9) sin(0)) (5) 

c = £(PC(0)sm(0)) 2 , (6) 



where PC{9) refers to the phase congruency value determined at orientation 0, 
and the sum is performed over the discrete set of orientations used (typically 
six). The angle of the principal axis # is given by 

The maximum and minimum moments, M and m respectively, are given by 

M=l{c + a+^V + {a-cY) (8) 
m = ±(c + a-^/b 2 + (a-c) 2 ). (9) 



This calculation of the maximum and minimum moments, along with the prin- 
cipal axis, corresponds to preforming a singular value decomposition on a phase 
congruency covariance matrix. The moments correspond to the singular values. 
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3.1 Comparison with the Harris Operator 

A comparison with the Harris operator is appropriate given its widespread 
use. The analysis described above is similar to that adopted by Harris and 
Stephens [3]. However they consider the minimum and maximum eigenvalues, 
a and /?, of the image gradient covariance matrix in developing their corner 
detector. The gradient covariance matrix is given by 

(10) 

where I x and I y denote the image gradients in the x and y directions. A 'corner' 
is said to occur when the two eigenvalues are large and similar in magnitude. To 
avoid an explicit eigenvalue decomposition Harris and Stephens devise a measure 
using the determinant and trace of the gradient covariance matrix 

R = det(G) - k(tr(G)) 2 , (11) 

where det(O) = and tr(G) = a + /?, the parameter k is traditionally set to 
0.04. This produces a measure that is large when both a and /? are large. However 
we have the problem of determining what is large. Noting that elements of the 
image gradient covariance matrix have units of intensity gradient squared we can 
see that the determinant, and hence the measure R will have units of intensity 
gradient to the fourth. This explains why the Harris operator is highly sensitive 
to image contrast variations which, in turn, makes the setting of thresholds 
exceedingly difficult. Some kind of sensitivity to image contrast is common to all 
corner operators that are based on the local autocorrelation of image intensity 
values and/or image gradient values. 

Unlike image intensity gradient values phase congruency values are normal- 
ized quantities that have no units associated with them. If the moments are 
normalized for the number of orientations considered we end up with phase con- 
gruency moment values that range between 0 and 1. Being moments these values 
correspond to phase congruency squared. Accordingly we can use the maximum 
and minimum phase congruency moments directly to establish whether we have a 
significant edge and/or corner point. It should be emphasized that the minimum 
and maximum moments provide feature information in their own right; one does 
not have to look at their ratios. We can define a priori what a significant value 
of phase congruency moment is, and this value is independent of image contrast. 



4 Results 



The performance of the phase congruency operator was compared to the Har- 
ris operator on a synthetic test image, and on a real scene containing strong 



shadows. The results are shown in figures 3 and 4. Raw Harris corner strength 
and raw phase congruency corner and edge strength images are displayed for 
comparison. It should be noted that the Harris corner strength values varied by 
many orders of magnitude across the image depending on contrast. To facilitate 
the display of the Harris corner strength image what is actually shown here is 
the fourth root of the image. Even after this transformation large variations are 
still evident. In contrast the phase congruency edge and corner strength images 
are minimally affected by image contrast and are readily thresholded (in this 
case with a value of 0.4 on both images) to produce a clear set of features. 

In applying the Harris operator to the synthetic image the standard deviation 
of the smoothing Gaussian was deliberately set to a large value of 3 pixels to 
illustrate the problem the operator has in localizing 'T' junctions. Notice how 
the detected Harris corners are displaced inwards on the left hand grey scale, 
note also the double response where the line intersects the circle. The phase 
congruency operator, on the other hand, locates T' junctions precisely. Changing 
the number of filter scales used to compute phase congruency does not affect the 
localization, it only affects the relative significance of features at different scales. 
For the real scene the standard deviation of the smoothing Gaussian for the 
Harris operator was 1 pixel. For both images phase congruency was computed 
using Gabor filters over 4 scales (wavelengths of 4, 8, 16 and 32 pixels) and over 
6 orientations. 

Another important point to note is that the phase congruency edge map 
includes the corner map, this is unusual for an edge detector! The fact that the 
phase congruency corner map is a strict subset of the phase congruency edge 
map greatly simplifies the integration of data computed from edge and corner 
information. This facilitates the process of building a model of the scene from 
point and edge data matched over two or more views. In contrast, if the edge and 
corner information is computed via separate means, say the Canny and Harris 
operators respectively, the edge data is unlikely to include the corner data. The 
Gaussian smoothing that is applied to reduce the influence of noise results in 
the Canny edges being weakened in strength, and rounded at corner locations. 

MATLAB code is available for those wishing to replicate the results presented 
here [21]. 

5 Conclusion. 

Phase congruency provides a contrast invariant way of identifying features within 
images. By combining phase congruency information over multiple orientations 
into a covariance matrix, and calculating the minimum and maximum moments 
we produce a highly localized operator that can be used to identify both edges 
and corners in a contrast invariant way. The contrast invariance facilitates the 
tracking of features over extended image sequences under varying lighting con- 
ditions. An additional advantage of the operator is that the phase congruency 
corner map is a strict subset of the phase congruency edge map. This simplifies 
the integration of data computed from edge and corner information. 



Original image. 



Fourth root of Harris corner 
strength image (<r = 3). 




Phase congruency edge strength 




Harris corners with threshold 10 7 
(maximum corner strength was 
3.4 x 10 9 ). 



O J 



Phase congruency corner strength 
image. 




Phase congruency corners with 
threshold 0.4 (maximum possible 
phase congruency value is 1). 



Fig. 3. Comparison of the Harris and phase congruency operators on a test image. 



Original image. Fourth root of Harris corner strength im- 

age (a = 1). 




Phase congruency edge strength image. Phase congruency corner strength image. 




Harris corners with threshold 10 s (maxi- Phase congruency corners with threshold 
mum corner strength was 1.25 x 10 10 ). 0.4 (maximum possible phase congruency 
value is 1). 



Fig. 4. Comparison of the Harris and phase congruency operators on an image with 
strong shadows. 
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Morphological Corner Detection 

This operator uses a variant of the morphological closing operator, which we have called asymmetrical closing. It consists of 
the successive application of different morphological tranformations using different structuring elements. Each of these 
structuring elements used to probe the image under study is tuned to affect corners of different orientation and brightness. 
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The operator can be easily implemented with Intel OpenCV as follows: 

cvDilate (src, tmpl, cross) ,- 
cvErode(tmpl,tmp2, lozeng) ,- 
cvDilate (src ,dst, ex) ,- 
cvErode (dst , tmpl , square) ; 
cvAbsDiff (tmp2, tmpl, dst) ; 
if (thresholds) 

cvThreshold(dst,det, threshold, 255, CV_THRESH_BINARY} ; 

Publications: 

Robert Laganiere, 
Morphological Corner Detection , 

in Proc. International Conference in Computer Vision, pp. 280-285, Bombay, India, January 1998. 
PDF [0.223 mb] 

Robert Laganiere, 

A Morphological Operator for Corner Detection , 

in Pattern Recognition, pp. 1643-1652, November 1998. 
Online version 

A modified version of this operator has been proposed in: 

Shih, Frank Y.; Chuang, Chao-Fa; Gaddipati, Vijayalakshmi, 
A modified regulated morphological corner detector , 
in Pattern Recognition Letters, pp. 931-937, April 2005. 
Online version 

Source code 
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Available here. The file comorph.h contains the operator's actual implementation that can be reused in other applications. 



JUDOCA: junction Detection Operator based on Circumferential 
Anchors 

Participant: Rimon Elias 

This junction detection operator defines junctions as points where linear ridges in the gradient domain intersect. The radial 
lines that compose the junction are therefore identified by searching, in a circular neighborhood, for directional maxima of the 
intensity gradient. The proposed algorithm operates on two binary edge maps, the computational complexity of the detection 
process is then considerably reduced. 

Algorithm Outline: 

1 . Apply vertical and horizontal Gaussian derivative filters on image I. 

2. Compute gradient magnitude and create two binary images from it: 

1 . The first one, B, is created from the imposition of a threshold, t B , on the gradient image. 

2. The second one, B + , contains the points of B that are local maxima in the direction of the gradient. 

3. For each point p in B, consider a circle of radius rho centered on this point and obtain the list of candidate points q t in 
B + that lie on the circumference of this circle (see Figure 1). These so-called circumferential anchor (CA) points are 
the extremities of potential radial lines for the putative junction. 




(a) <b) 
Figure 1 : (a) One corner of a box that produces a Y-junction. (b) The junction at p with three circumferential anchors 
q p q 2 and q 3 (superimposed on the gradient image). 

4. For each CA point q t in the list, consider the set of all points located at a distance iess than one pixel to the segment 
that joins the current CA point to the central point p. This set is used to determine if the corresponding putative 
junction radial line will be accepted (considering the points in B) and if yes, what the strength of this junction line would 
be (considering the points in B + ), that is: 

1 . To be accepted as a junction line, a continuous path of B points joining the CA point and the central point p 
must exist. If not, then reject this radial line and repeat the scanning operation with the next CA point. 

2. The strength, SJ<p,q>), of this junction radial line is defined as the sum of the squared distances from the B + 
points in the currently considered set to the <p,q> line segment. This strength is normalized by the length of 

~ this segment. 

5. If the strength of this junction radial line is smaller than a predetermined threshold, s Jt then reject this radial line. 
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Otherwise, <p,q> becomes one of the branch of the putative junction at p. 
6. If the number of branches found at p is less than 2, then there is no junction at this location. Otherwise, record the 
orientation of the accepted junction radial lines and set the strength, S/p), of this junction as being the minimum 
across all radial line strengths, i.e.: 

SJ P )=miniS/< P ,q>) 



Examples: 




(Click on image for better resolution) 




(Click on image for better resolution) 



Publications: 

Robert Laganiere, Rimon Elias 

The Detection of Junction Features in Images , 

in Proc. international Conference on Acoustic, Speech and Signal processing, pp. 573-576, Montreal, Canada, May 2004. 
PDF [0.223 mb] 

Source code 

Available here . 



Copyright © 2004 VIVA Lab 
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Corner Detection 



Functions 

template<...> void cornerResponseFunction (Srclterator sul, Srclterator sir, SrcAccessor as, 
Destlterator dul, DestAccessor ad, double scale) 

Find corners in an image (1). 

template<...> void foerstnerCornerDetector (Srclterator sul, Srclterator sir, SrcAccessor as, 
Destlterator dul, DestAccessor ad, double scale) 

Find corners in an image (2). 

template<...> void rohrCornerDetector (Srclterator sul, Srclterator sir, SrcAccessor as, Destlterator dul, 
DestAccessor ad, double scale) 

Find corners in an image (3). 

template<...> void beaudetCornerDetector (Srclterator sul, Srclterator sir, SrcAccessor as, Destlterator 
dul, DestAccessor ad, double scale) 

Find comers in an image (4). 



Detailed Description 

Measure the 'cornerness' at each pixel. Note: The Kitchen-Rosenfeld detector is not implemented 
because of its inferior performance. The SUSAN detector is missing because it's patented. 



Function Documentation 

void beaudetCornerDetector (...) 

Find corners in an image (4). 

This algorithm implements a corner detector according to [P.R. Beaudet: "Rotationally Invariant 
Image Operators", Proc. Intl. Joint Conf. on Pattern Recognition, Kyoto, Japan, 1978, pp. 579-583]. 

The algorithm calculates the corner strength as the negative determinant of the Hessian Matrix. The 
local maxima of the corner strength denote the corners in the gray level image. 

The source value type must be a linear algebra, i.e. addition, subtraction, and multiplication with itself, 
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muitiplication with doubles and NumericTraits must be defined. 
Declarations: 

pass arguments explicitly: 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

void 

beaudetCornerDetector (Srclterator sul, Srclterator sir, SrcAccessor as, 
Destlterator dul, DestAccessor ad, 
double scale) 

} 

use argument objects in conjunction with Argument Object Fa ctories: 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

inline 

void beaudetCornerDetector ( 

triple<Srclterator, Srclterator, SrcAccessor> src, 
pair<DestIterator, DestAccessor> dest, 
double scale) 

} 

Usage: 

#inciude " vigra/cornerdetection.hxx " 
Namespace: vigra 

vigra : : BImage src ( w, h) , corners (w, h) ; 
vigra: :FImage beaudet_corner_strength Cw,h) ; 

// empty corner image 
corners . init (0.0) ; 



// find corner response at scale 1.0 

vigra: : beaudetCornerDetector (srcImageRange (src) , destlmage (beaudet_corner_strength) , 
1.0); 

// find local maxima of corner response, mark with 1 

vigra: : localMaxima (srcImageRange (beaudet_corner_strength) , destlmage (corners) ) ; 

Required Interface: 

Srclmagelterator src__upperlef t, src_lowerright; 
Destlmagelterator dest_upperlef t ; 

SrcAccessor src_accessor ; 
DestAccessor dest_accessor; 

SrcAccessor: :value_type u = src_accessor (src_upperlef t) ; 
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double d; 

u = u + u 
u = u - u 
u = u * u 
u = d * u 

dest_accessor. set (u, dest_upperlef t) ; 



void cornerResponseFunction (...) 



Find corners in an image (1 ). 

This algorithm implements the so called 'corner response function' to measure the 'cornerness' of 
each pixel in the image, according to [C.G. Harris and M.J. Stevens: "A Combined Corner and Edge 
Detector", Proc. of 4th Alvey Vision Conference, 1988]. Several studies have found this to be a very 
robust corner detector, although it moves the corners somewhat into one region, depending on the 
scale. 

The algorithm first determines the structure tensor at each pixel by calling structureTensor(). Then 
the entries of the structure tensor are combined as 

CoruerltespoHSC « det{StryctuwTensor)-a04tr(StructurcsT«ufs>r) 3 = AB-C^-OMiA+B) 2 

The local maxima of the corner response denote the corners in the gray level image. 

The source value type must be a linaer algebra, i.e. addition, subtraction, and multiplication with itself, 
multiplication with doubles and Nume ri cTraits must be defined. 

Declarations: 

pass arguments explicitly: 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

void 

cornerResponseFunction (Srclterator sul, Srclterator sir, SrcAccessor as, 
Destlterator dul, DestAccessor ad, 
double scale) 

} 

use argument objects in conjunction with Argumen t Object Factories : 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

inline 

void cornerResponseFunction ( 

triple<SrcIterator, Srclterator, SrcAccessor> src, 
pair<DestIterator, DestAccessor> dest, 
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double scale) 

} 

Usage: 

#include "vigra/ cornerdetection.hxx " 
Namespace: vigra 

vigra: :BImage src(w,h), corners (w,h) ; 
vigra: :FImage corner_response (w, h) ; 

// empty corner image 
corners . init (0.0) ; 



// find corner response at scale 1.0 

vigra: : cornerResponseFunction (srcImageRange (src) , destlmage (corner_response) , 
1.0); 

// find local maxima of corner response, mark with 1 

vigra : : localMaxima (srcImageRange (corner_response) , destlmage {corners) ) ; 

// threshold corner response to keep only strong corners (above 400.0) 
transformlmage (srcImageRange (corner_response) , destlmage (corner_response) , 
vigra: :Threshold<double, double >{ 

400.0, std: :numeric_limits<double> : :max() , 0.0, 1.0)); 

// combine thresholding and local maxima 

vigra: : comb in eTwo Images (srcImageRange (corners) , srclmage (corner_response) , 
destlmage (corners) , std: :multiplies<f loat> () ) ; 

Required Interface: 

Srclmagelterator src_upperlef t, src_lowerright; 
Destlmagelterator dest__upperlef t ; 

SrcAccessor src_accessor ; 
DestAccessor dest_accessor 

SrcAccessor: :value_type u = src_accessor (src_upperlef t) ; 
double d; 

u = u + u 
u = u - u 
u = u * u 
u = d * u 

dest_accessor.set (u, dest_upperlef t) ; 



void foerstnerCornerDetector (...) 



Find corners in an image (2). 

This algorithm implements the so called 'Foerstner Corner Detector" to measure the 'cornerness' of eac 
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pixel in the image, according to [W. Forstner: "A feature based correspondence algorithms for image 
matching", Intl. Arch. Photogrammetry and Remote Sensing, vol. 24, pp 160-166, 1986]. It is also know 
as the "Plessey Detector" by Harris. However, it should not be confused with the "Corner Repsonse 
Function ", another detector invented by Harris. 

The algorithm first determines the structure tensor at each pixel by calling structureTensorf). Then th< 
entries of the structure tensor are combined as 

„ det(StractureTeasor) AB-G 2 

FoerstnerCornerStrength = tr(StmctureTeu&or) = -X+W 

The local maxima of the corner strength denote the corners in the gray level image. Its performance is 
similar to the cornerResponseFunction(). 

The source value type must be a division algebra, i.e. addition, subtraction, multiplication, and division 
with itself, multiplication with doubles and Nu mericTraits must be defined. 



Declarations: 



pass arguments explicitly: 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

void 

foerstnerCornerDetector {Srclterator sul, Srclterator sir, SrcAccessor as, 
Destlterator dul, DestAccessor ad, 
double scale) 

} 



use argument objects in conjunction with Argument O bject Factories: 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

inline 

void foerstnerCornerDetector ( 

triple<SrcIterator, Srclterator, SrcAccessor> src, 
pair<DestIterator, DestAccessor> dest, 
double scale) 



Usage: 

#include "vigra/c ornerdetection.hxx " 
Namespace: vigra 

vigra: :BImage src (w,h) , corners (w,h) ; 

vigra : : Flmage f oerstner_corner_strength (w, h) ; 



If empty corner image 
corners. init (0.0) ,- 
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// find corner response at scale 1.0 

vigra: : f oerstnerCornerDetector (srcImageRange (src) , destlmage (f oerstner_corner_strength) 
1.0); 

// find local maxima of corner response, mark with 1 

vigra: : localMaxima (srcImageRange (foerstner_corner_strength) , destlmage (corners) ) ; 

Required Interface: 

Srclmagelterator src_upperlef t , src_lowerright; 
Destlmagelterator dest_upperlef t ; 

SrcAccessor src_accessor ; 
DestAccessor dest_accessor; 

SrcAccessor: :value_type u = src_accessor (src_upperlef t) ; 
double d; 

u = u + u 
u = u - u 
u = u * u 
u = u / u 
u = d * u 

dest_accessor.set (u, dest_upperlef t) ; 



void rohrCornerDetector (...) 



Find comers in an image (3). 

This algorithm implements yet another structure tensor-based corner detector, according to [K. Rohr: 
"Untersuchung von grauwertabhangigen Transformationen zur Ermittlung der optischen Flusses in 
Bildfolgen", Diploma thesis, Inst, fur Nachrichtensysteme, Univ. Karlsruhe, 1987, see also K. Rohr: 
"Modelling and Identification of Characteristic Intensity Variations", Image and Vision Computing 10:2 
(1992) 66-76 and K. Rohr: "Localization Properties of Direct Corner Detectors", J. of Mathematical 
Imaging and Vision 4:2 (1994) 139-150]. 

The algorithm first determines the structure tensor at each pixel by calling structureTensor(). Then 
the entries of the structure tensor are combined as 

RobrCoraetStrengtli = det(Struc*urvTtm$<>r) = AB — C a 

The local maxima of the corner strength denote the corners in the gray level image. Its performance 
is similar to the cornerResponseFunction(). 

The source value type must be a linear algebra, i.e. addition, subtraction, and multiplication with itself, 
multiplication with doubles and NumericTraits must be defined. 

Declarations: 
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pass arguments explicitly: 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

void 

rohrCornerDetector (Srclterator sul, Srclterator sir, SrcAccessor as, 
Destlterator dul, DestAccessor ad, 
double scale) 

} 

use argument objects in conjunction with Argument Object Factories : 

namespace vigra { 

template <class Srclterator, class SrcAccessor, 

class Destlterator, class DestAccessor> 

inline 

void rohrCornerDetector ( 

triple<SrcIterator, Srclterator, SrcAccessor> src, 
pair<DestIterator, DestAccessor> dest, 
double scale) 



Usage; 

#include "vig ra/cornerd etec tion.hxx " 
Namespace: vigra 

vigra: :BImage src(w,h), corners (w, h) ; 
vigra: :FImage rohr_corner_strength{w, h) ; 

// empty corner image 
corners . init (0.0) ; 



// find corner response at scale 1.0 

vigra: :rohrCornerDetector(src ImageRange (src), de s t Image ( r ohr_c orne r_s t r ength ) , 
1.0); 

// find local maxima of corner response, mark with 1 

vigra : : localMaxima (srdmageRange (rohr_corner_strength) , destlmage (corners) ) ; 

Required Interface: 

Srclmagelterator src_upperlef t, src_lowerright; 
Destlmagelterator dest_upperlef t ; 

SrcAccessor src_accessor ; 
DestAccessor dest_accessor; 

SrcAccessor: :value_type u = src_accessor (src_upperlef t) ; 
double d; 

u = u + u 
u = u - u 
u = u * u 
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u = d * u 

dest_accessor. set (u, dest_upperlef t) ; 
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Next: Line segment detection Up: Feature detection Previous: Dis playing a corner map Contents . 

The Harris corner detector 

#include <gandalf /vision/harris_corner.h> 

The Harris corner detector [#! Harris: Stephens :ALVEY88!#] computes the locally averaged moment matrix 
computed from the image gradients, and then combines the eigenvalues of the moment matrix to compute a corner 
"strength", of which maximum values indicate the comer positions. Here is an example code fragment using the 
Harris corner detector. 

Gan_lmage *plmage; /* declare image from which corners will be computed */ 
Gan_MasklD *pFilter; /* convolution mask */ 
Gan_CornerFeatureMap CornerMap; /* declare corner map */ 

/* . . . fill image ... */ 

/* initialise corner map */ 
gan_corner_feature_map_form ( fcCornerMap, 

1000 ) ; /* initial limit on number of corners */ 

/* create convolution mask */ 

pFilter = gan_gauss_mask_new ( GAN_FL0AT, 1.0, 9, 1.0, NULL ); 

/* apply Harris corner detector */ 

gan_harris_corner_g ( plmage, /* input image */ 

NULL, /* or binary mask of pixels to be processed */ 

NULL, NULL, /* or image pre -smoothing masks */ 

pFilter, pFilter, /* gradient smoothing */ 

0.04, /* kappa used in computing corner strength */ 

0.04, /* corner strength threshold */ 

NULL, /* or affine coordinate transformation */ 

0, /* status value to assign to each corner */ 

NULL, /* or pointer to camera structure defining 

distortion model */ 
NULL, /* or parameters of local feature map */ 
&CornerMap ) ; /* result corner map */ 

/* free convolution mask and corner map */ 

gan_masklD_free ( pFilter ) ; 

gan_corner_f eature_map_f ree { fcCornerMap ) ; 
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Abstract. A critical step in automatic fingerprint recognition is the accurate seg- 
mentation of fingerprint images. The objective of fingerprint segmentation is to 
decide which part of the images belongs to the foreground containing features 
for recognition and identification, and which part to the background with the 
noisy area around the boundary of the image. Unsupervised algorithms extract 
blockwise features. Supervised method usually first extracts point features like 
coherence, average gray level, variance and Gabor response, then a Fisher linear 
classifier is chosen for classification. This method provides accurate results, but 
its computational complexity is higher than most of unsupervised methods. This 
paper proposes using Harris corner point features to discriminate foreground and 
background. Shifting a window in any direction around the corner should give a 
large change in intensity. We observed that the strength of Harris point in the fore- 
ground area is much higher than that ofHarris point in background area. The un- 
derlying mechanism for this segmentation method is that boundary ridge endings 
are inherently stronger Harris corner points. Some Harris points in noisy blobs 
might have higher strength, but it can be filtered as outliers using corresponding 
Gabor response. The experimental results proved the efficiency and accuracy of 
new method are markedly higher than those of previously described methods. 

1 Introduction 

The accurate segmentation of fingerprint images is key component to achieve high per- 
formance in automatic fingerprint recognition systems. If more background areas are 
included into segmented fingerprint of interest, more false features are possibly intro- 
duced into detected feature set; If some parts of foreground are excluded, useful feature 
points may be missed. There are two types of fingerprint segmentation algorithms: un- 
supervised and supervised. Unsupervised algorithms extract blockwise features such as 
local histogram of ridge orientation [1,2], gray-level variance, magnitude of the gra- 
dient in each image block [3], Gabor feature [4,5]. Practically, the presence of noise, 
low contrast area, and inconsistent contact of a fingertip with the sensor may result in 
loss of minutiae or more spurious minutiae. Supervised method usually first extracts 
several features like coherence, average gray level, variance and Gabor response [5-7], 
then a simple linear classifier is chosen for classification. This method provides accurate 
results, but its computational complexity is higher than most unsupervised methods. 

Segmentation in low quality images faces several challenging technical problems. 
First problem is the presence of noise caused by dust and grease on the surface of live- 
scan fingerprint scanners. Second problem is ghost images of fingerprints remaining 



from the previous image acquisition [7]. Third problem is low contrast fingerprint ridges 
generated through inconsistent contact press or dry/wet finger surface. Fourth problem 
is indistinct boundary if the features in the fixed size of window are used. Final problem 
is segmentation features being sensitive to the quality of image. 

This paper proposes using Harris corner point features [8, 9] to discriminate fore- 
ground and background. The Harris corner detector was developed originally as fea- 
tures for motion tracking, it can reduce significantly amount of computation compared 
to tracking every pixel. It is translation and rotation invariant but not scale invariant. We 
found that the strength of a Harris point in the foreground area is much higher than that 
of a Harris point in the background area. Some Harris points in noisy blobs might have 
higher strength, but it can be filtered as outliers using corresponding Gabor response. 
The experimental results proved the efficiency and accuracy of new method are much 
better than those of previously described methods. Furthermore, this segmentation al- 
gorithm can detect accurate boundary of fingerprint ridge regions, which is very useful 
in removing spurious boundary minutiae, and most current segmentation methods can 
not provide consistent boundary minutiae filtering. 

2 Features for Fingerprint Segmentation 

Feature selection is the first step for designing fingerprint segmentation algorithm. 
There are two general types of features used for fingerprint segmentation, i.e., block 
features and pointwise features. In [5, 10] selected point features include local mean, 
local variance or standard deviation, and Gabor response of the fingerprint image. Lo- 
cal mean is calculated as Mean = J^w ^ ^ oca ^ var i ance is calculated as Var = 
Ylwi 1 ~ Mean) 2 , w is window size centered the processed pixel. The Gabor response 
is the smoothed sum of Gabor energies for eight Gabor filter responses. Usually the 
Gabor response is higher in the foreground region than that in the background region. 
The coherence feature indicates how strong the local window gradients centered the 
processed point in the same dominant orientation. Usually the coherence will be much 
higher in the foreground than in the background, but it may be influenced significantly 
by boundary signal and noisy signal. Therefore, single coherence feature is not suffi- 
cient for robust segmentation. Systematic combination of those features is necessary. 



Because pointwise-based segmentation method is time consuming, blockwise fea- 
tures are usually used in the commercial automatic fingerprint recognition systems. 
Block mean, block standard deviation, block gradient histogram [1,2], block average 
magnitude of the gradient [11] are most common block features for fingerprint seg- 
mentation. In [12] gray-level pixel intensity-derived feature called block clusters de- 
gree(CluD) is introduced. CluD measures how well the ridge pixels are clustering. 




0) 



E{x,y) = \og{JJ g \F{r,d)\ 2 } 



(2) 



Texture features, such as Fourier spectrum energy [6], Gabor features [4, 1 3] and Gaussian- 
Hermite Moments [14], have been applied to fingerprint segmentation. Ridges and val- 
leys in a fingerprint image are generally observed to possess a sinusoidal-shaped plane 
wave with a well-defined frequency and orientation [15], and non-ridge regions does 
not hold this surface wave model. In the areas of background and noisy regions, it is 
assumed that there is very little structure and hence very little energy content in the 
Fourier spectrum. Each value of energy image E(x,y) indicates the energy content of 
the corresponding block. The fingerprint region may be differentiated from the back- 
ground by thresholding the energy image. The logarithm values of the energy is used 
to convert the large dynamic range to a linear scale(Equation 2). The region mask is 
obtained by thresholding E(x, y). However, uncleaned trace finger ridges and straight 
stripes are unfortunately included into regions of interest (Figure 1(c)). 




(a) (b) 




(c) (d) 




(e) (0 

Fig. 1. (a) and (b) arc two original images, (c) and (d) are FFT energy maps for images (a) and 
(b), (e) and (f) are Gabor energy maps for images (a) and (b), respectively 



Gabor filter-based segmentation algorithm is now most often used method [4, 13]. 
An even symmetric Gabor filter has the following spatial form: 



9(x, V, 0, /, a s , a y ) = ex V {-\\% + ^ftcos^fx,) (3) 

For each block of size W x W centered at (x,y), 8 directional Gabor features are com- 
puted for each block, the standard deviation of 8 Gabor features is utilized for segmen- 
tation. The formula for calculating the magnitude of Gabor feature is defined as, 

|(W2)-i («/a)-i 

G(X,Y,e,f,a.,a v )'-\ £ Yl HX + x 0 ,Y + yo)g(xo,yo,O l f,<r tt ,tr y ) 




However, fingerprint images with low contrast or false traces ridges or noisy com- 
plex background can not be segmented correctly by Gabor filter-based method (Figure 
1(e)). 

In [14], similarity is found between Hermite moments and Gabor filter. Gaussian- 
Hermite Moments has been successful ly used to segment fingerprint images in [ 14] . Or- 
thogonal moments use orthogonal polynomials as transform kernels and produce min- 
imal information redundancy, Gaussian Hermite moments(GHM) can represent local 
texture feature without minimal noise effect. 



3 Harris Corner Points 



3.1 Review of Harris-corner-points 

We propose using Harris corner point features [8, 9] to discriminate foreground and 
background. The Harris corner detector was developed originally as features for motion 
tracking, it can reduce significantly amount of computation compared to tracking every 
pixel. Shifting a window in any direction around the corner should give a large change 
in intensity. Comer points provide repeatable points for matching, so some efficient 
methods have been designed [8, 9]. Gradient is ill defined at a corner, so edge detectors 
perform poorly at comers. However, in the region around a comer, gradient has two 
or more different values. The comer point can be easily recognized by looking a small 
window. Shifting a window around a comer point in any direction should give a large 
change in gray-level intensity, 

Given a point I(x,y), and a sh.iR(Ax, Ay), the auto-correlation function E is defined 

as: 

e(x, y )= Yl [ j o*> w) - ^ + ^ » + A y^ 2 ® 

w(x,y) 

where w(x,y)is window function centered on image point(x,y). For a small sh\fis,[Ax,Ay\ , 
the shifted image is approximated by a Taylor expansion truncated to the first order 
terms, 

I(xi + Ax, Vi + Ay) t* I{xi,yi) + [I x (xi,yi)I v (xi,yi)] (6) 



where I x (xi,yi) and I y {x u y t ) denote the partial derivatives in x and y, respectively. 
Substituting approximation Equation 6 into Equation 5 yields, 

E{x,y)~Y, W {x,y)[i(xi>yi)- I ( x i+ As: >yi+ A y)] 2 

= Z Mx , y) (l(xi,Vi) ~ H^uVi) ~ [I*(*i,Vi)Iv(*i>Vi)\ [X]) 
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= [Ax,Ay]M(x,y) 



E{Ax,Ay) = [Ax,Ay]M(x,y)^ 



where M(x,y) is a 2 x 2 matrix computed from image derivatives, called auto-correlation 
matrix which captures the intensity structure of the local neighborhood. 



3 .2 Strength of Harris-corner points of a Fingerprint Image 

In order to detect interest points, the original measure of corner response in [8] is : 

tfet(M) AiA 2 fl0 > 
Trace(M) Ai + A 2 

The auto-correlation matrix (M) captures the structure of the local neighborhood. Based 
on eigenvalues^!, A 2 ) of M, interest points are located where there are two strong eigen 
values and the corner strength is a local maximum in a 3 x 3 neighborhood. To avoid the 
explicit eigenvalue decomposition of M, Trace{M) is calculated as 2| + 1%, Det(m) 
is calculated as - {I x I y ) 2 , and 

R = Det(m) -kx Trace(M) 2 (11) 

To segment the fingerprint area (foreground) from the background, the following 
"corner strength" measure is used, because there is one undecided parameter k in equa- 
tion(ll). 



3.3 Harris-corner-points based Fingerprint Image Segmentation 



We found that the strength of a Harris point in the fingerprint area is much higher than 
that of a Harris point in background area, because boundary ridge endings inherently 
possess higher corner strength. Most high quality fingerprint images can be easily seg- 
mented by choosing appropriate threshold value. In Figure 2, a corner strength of 300 
is selected to distinguish corner points in the foreground from those in the background. 
Convex hull algorithm is used to connect harris comer points located in the foreground 
boundary. 




(a) (b) (c) (d) (e) 

Fig. 2. a fingerprint with harris comer strength of (b)10, (c)60, (d)200, and (e)300. This fingerprint 
can be successfully segmented using corner response threshold of 300. 

It appears relatively easy for us to segment fingerprint images for following image 
enhancement, feature detection and matching. However, two technical problems need 
to be solved. First, different "corner strength" thresholds are necessary to achieve good 
segmentation results for different qualities images based on image characteristical anal- 
ysis. Second, Some Harris points in noisy blobs might have higher strength, it can not 
be segmented by choosing simply one threshold. When single threshold is applied to 
all the fingerprint images in one whole database, not all the corner points in the back- 
ground in a fingerprint image are removed, some corner points in noisy regions can not 
be thresholded even using high threshold value (Figure 3). In order to deal with such 
situations, we implemented a heuristic selection algorithm using corresponding Gabor 
response (Figure 4). 



4 Experimental Results 

The proposed methodology is tested on FVC2002 DB1 and DIM, each database consists 
of 800 fingerprint images (100 distinct fingers, 8 impressions each). Image size is 374 x 
388 and the resolution is 500dpi. To evaluate the methodology of adapting a gaussian 
kernel to the local ridge curvature of a fingerprint image, we modified Gabor-based 
fingerprint enhancement algorithm [15, 16] with two kernel sizes: the smaller one in 
high-curvature regions and the larger one in pseudo-parallel ridge regions, minutiae are 



(a) (b) (c) <d) 

Fig. 3. a fingerprint with harris corner strength of (a)100, (b)500, (c)1000, (d) 1500 and (e)3000. 
Some noisy corner points can not be filtered completely even using comer response threshold of 
3000. 



Fig. 4. Segmentation result and final feature detection result for the image shown in the Figure 
1(a). (a) Segmented fingerprint marked with boundary line, (b) final detected minutiae. 



detected using chaincode-based contour tracing [17], the fingerprint matcher developed 
by Jea et al. [18] is used for performance evaluation. 

Our methodology has been tested on low quality images from FVC2002 . To validate 
the efficiency of proposed segmentation method, current widely-used Gabor filter-based 
segmentation algorithm [4, 13] and NIST segmentation [19] are utilized for compari- 
son. 

The proposed segmentation method have a remarkable advantage over current meth- 
ods in terms of boundary spurious minutiae filtering. Figure 5 (a) and (b) show unsuc- 
cessful boundary minutiae filtering using NIST method [19], which is implemented by 
removing spurious minutiae pointing to invalid block and removing spurious minutiae 
near invalid blocks, and invalid blocks are defined as blocks with no detectable ridge 
flow. However, boundary blocks are more complicated, so the method in [19] fails 
to remove most boundary minutiae. In Figure 5 (c) and (d) show the filtering results 
of proposed method. In comparison of Figure 5(a) against (c) and (b) and (d), 30 and 
17 boundary minutiae are filtered, respectively. Performance evaluations for FVC2002 
DB1 and DB4 are shown in Figure 6. For DB1, ERR for false boundary minutiae filter- 
ing using proposed segmented mask is 0.0 1 06 and EER for NIST boundary Filtering is 
0.0125. For DB4, ERR for false boundary minutiae filtering using proposed segmented 
mask is 0.0453 and EER for NIST boundary Filtering is 0.0720. 




Fig. 6. ROC curves for (a) FVC2002 DB1 and (b) FVC2002 DB4. 



5 Conclusions 



In this paper, a robust interest point based fingerprint segmentation is proposed for 
fingerprints of varied image qualities. The experimental results compared with those 
of previous methods validate that our algorithm has better performance even for low 
quality images, in terms of including less background and excluding less foreground. In 
addition, this robust segmentation algorithm is capable of filtering efficiently spurious 
boundary minutiae. 
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Abstract 

Noise filtering of images is essentially a smoothing process, and it is an issue that has been addressed for many years. The most 
commonly used low-pass filtering methods blur important image structures such as edges and lines, and thus reduce image contrast 
and damage image fidelity. This paper presents a structure adaptive anisotropic filtering technique with its application to processing 
magnetic resonance images. It differs from other techniques in that, instead of using local gradients as a means of controlling the 
anisotropism of filters, it uses both a local intensity orientation and an anisotropic measure of level contours to control the shape and 
extent of the filter kernel. This ensures that corners and junctions are well preserved throughout the filtering process. The following 
two aspects of the proposed technique demonstrate the advantage of using this filtering method, Firstly, the use of local orientation 
detection provides a robust and convenient way for shaping the filter kernel. Secondly, the structural adaptiveness of the filtering 
process ensures that salient image features are non-symmetrically enhanced. 

Keywords: Adaptive filtering; Orientation detection; Magnetic resonance imaging 



1. Introduction 

In medical imaging one is often faced with the problem 
of compromising between image quality and data acqui- 
sition time. To achieve a good Signal-to-Noise Ratio 
(SNR), time domain averaging is usually used. This 
prolongs the imaging time, and the method itself is not 
suited for studies that rely on temporal resolution to 
reveal important anatomical and functional information. 
While providing an ultimate solution to the problem 
mentioned above, improvement of acquisition hardware 
is often restrained by the practicality and costs involved. 
In this situation, post processing of image data is needed 
to improve the clarity of object details. Over the years, 
two broad categories of technique, linear and non-linear 
filtering methods, have been used for this purpose. The 
linear approach has the advantage of simplicity, but it 
does not respect regional structure and the resulting 
images often appear blurry and diffused. This 
undesirable effect can be reduced to a certain extent by 
the use of the non-linear scheme in which local structures 
and statistics are taken into account during the filtering 
process [1-3]. Because of their adaptiveness to local 
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image features, they are referred to as the adaptive 
methods. 

The idea of adaptive smoothing itself is not new, and 
many different approaches have been proposed over the 
years. A detailed overview and evaluation of some of the 
earlier methods can be found in [4], Lev, Zucker and 
Rosenfeld [5] performed one of the earliest studies in 
this field. They suggested applying at each image point 
a weighted mask with coefficients that are based on the 
differences between the value at the centre point and 
those of its neighbours. Similar approaches can also be 
found elsewhere [2,6]. The major drawback of these 
iterative smoothing methods is that their convergence 
properties are difficult to establish. A more recent 
approach, proposed by Blake and Zisserman [7], used 
weak continuity constraints to allow discontinuities in 
a piecewise continuous reconstruction of a noisy signal. 
The convergence behaviour was well studied and results 
provided were promising, but unfortunately it was 
accomplished with a formidable computational cost. 
By casting the problem in terms of a heat equation in 
an anisotropic medium, Perona and Malik [8] presented 
an anisotropic diffusion scheme. The principal aim of 
their method is to enhance edges and facilitate region 
segmentation, however, when small iteration steps are 
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used it can also be considered as a useful tool for image 
noise filtering [9]. 

When developing a filtering method for image data, 
the following two basic requirements should generally be 
observed: efficiency and fidelity. Efficiency requires noise 
to be effectively removed and the morphological defini- 
tion of the image to be enhanced, whereas fidelity 
demands detailed image structures to be preserved and 
no artifacts are generated during the filtering process. 
Based on these two criteria, in this paper we introduce 
a structure adaptive anisotropic filtering scheme with its 
application to processing magnetic resonance images. It 
differs from other techniques in that, instead of using 
local gradients as a means of controlling the anisotrop- 
ism of the filters, it uses the local intensity orientation of 
level contours and an anisotropic measure to control the 
shape and the extent of the filter kernel. The relationship 
between the proposed technique and other methods will 
be discussed, and a qualitative evaluation of the algo- 
rithm will be performed by using a synthetic image 
with various degrees of additive noise. 



2. Basic principle of the algorithm 

The most commonly used filtering algorithm computes 
a new image f s (x) by applying, at each point * 0 , a kernel 
k(x 0 , x) to the original image f(x) so that: 

where: 

is a normalization factor. In Eq. (1) x = (JC|,»2, . . . , x m ) 
is an m dimensional position vector. It is worth mention- 
ing at this stage that the basic principle of the algorithm 
proposed is applicable to m dimensional pictures, but 
applications provided in this paper will be restricted to 
2D and 3D magnetic resonance images. When function k 
in Eq. (1) is invariant over space, it uniformly smooths 
out image noise without paying any attention to image 
details. The basic idea behind the proposed algorithm is 
that the kernel k{x 0 ,x) used in Eq. (1) should be made 
variable and allow to be shaped or scaled according to 
local features within the neighbourhood U of x 0 . If one 
takes local gradients as the primary concern, then 
intuitively one requires the kernel to be broad where 
the gradient is weak so that uniform regions can be 
smoothed out. On the other hand, at piaces where the 
gradient is strong the kernel should be narrow so that 
edges and corners are preserved. 

When a Gaussian kernel is adopted, the relationship 
stated above can be expressed analytically by the 



following equation: 

k(x 0 , x) = p(x - * 0 ) e-l*-*l 2 /<*l v '<*>l> (3) 
where p(x) is a positive and rotationally symmetric cut- 
off function that satisfies the condition p(x) = 1 when 
jjrj < r, and r is the maximum support radius. The effec- 
tive support radius of the kernel is actually controlled by 
cr(x), which is a decreasing function with respect to 
|V/(* 0 )|- the local gradient strength. It is not difficult 
to see the drawback associated with this approach. 
Because the effective support radius is governed only 
by the local gradient strength, the filtering process is 
practically brought to a halt when an edge is encoun- 
tered. A better approach would be to shape the kernel 
so that it is only narrow in the predominant gradient 
direction and a reasonable effective support radius is 
maintained to smooth along but not across the edge. In 
this case, the shape of the kernel becomes elliptical as 
shown in Fig. 1, where the axis that lies along the edge 
is called the principal (major) axis. 

To accommodate this anisotropic nature of the filter- 
ing process, the following modification to Eq. (3) is 
necessary: 

k(x 0) x) = p{x - ^e-W'-^-^^^WC-^'-"^^)} 

(4) 

where n and n L are mutually normal unit vectors, and n is 
in parallel with the principal axis. In Eq. (4), (x -x Q )-n 
and (jc - jr 0 ) • n L are the dot products between vectors 
(x - jc 0 ) and n or n x , respectively. The shape of the kernel 
is controlled through two non-negative functions er^x) 
and <r 2 (x). A relationship of o\{x) +<rl{x) = ^(x) and 
a\{x) ~£ a 2 (x) is usually maintained so that k(x 0 ,x) is 




Fig. 1. Shaping of the filter kernel. When an edge is encountered the 
kernel is deformed into an ellipse with its principal axis aligned in 
parallel with the edge, and the smoothing is performed along but not 
across the edge. 



c 
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elongated along n. One may notice that compared with Eq. 
(3), the constraints imposed on o\{x) and a 2 {x) in Eq. (4) 
have been loosened. They are no longer governed directly 
by the local gradient strength and, as we shall see later, 
they are also influenced by corners and other local image 



Fig. 2. Several examples of orientated patterns (a)-(c) and their corresponding Fourier transform <e)-(g). They demonstrate that the power spectrum 
of such patterns clusters along a line through the origin in Fourier domain. The direction of the line is perpendicular to the dominant spatial 
n-orientated pattern (d), the power spectrum (h) spreads out in all directions. 

Bigun and Granlund [10] and that of Kass and Witkin 
[1 1], The level contours (isogrey values) of an orientated 
pattern consist of parallel lines. This definition for an 
orientated pattern includes not only different forms of 
edges (step, slope and second order edges, etc.), but 
also patterns, as shown in Fig. 2. The power spectrum 
of such a pattern clusters along a line through the origin 
in the Fourier domain, and the direction of the line is 
perpendicular to the dominant spatial orientation. If we 
relate the infinitesimal energy \F(u)\ 2 du of the Fourier 
transform to a mass distribution, then the second 
moment (the moment of inertia of \F(u)\ 2 du), E(n p ), 
with respect to a line through the origin with a unit 
normal vector n p is: 



To use Eq. (4) we have to define the direction of 
vectors n and n ± . The solution appears to be simple at 
a first thought, as one can define n to be perpendicular to 
the local gradient direction. And, indeed, it was the 
approach that was adopted initially. Problems occur 
when the signal-to-noise ratio of the image is low and 
its derived local gradient vectors are noisy. To circum- 
vent this problem, one can convolve the image with a 2D 
Gaussian kernel prior to calculating local intensity gra- 
dients. However, the smaller the signal-to-noise ratio, the 
larger the convolution window should be used. In this 
case, one is unintentionally assuming that the original 
signal varies smoothly within the selected window, and 
this assumption inevitably introduces errors when fine 
image details are encountered. If the orientation of the 
kernel is not properly aligned for anisotropic filtering, 
the image produced could, in fact, be worse than that 
of using a simple linear filtering technique. This 
motivated us to search for a better way of defining the 
filter kernel. 



3. Local orientation detection 



The basic idea behind our chosen technique for calcu- 
lating the direction of n was inspired by the work of 



E{n p )=\ i ...\d 2 p {u)\F{u)\ 2 du 



(5) 



where d p (u) is a real valued function which gives the 
shortest Euclidean distance between a point u and a 
line that is perpendicular to n p , i.e.: 

d 2 p {u) = (« - n p ) 2 = (nJu)(u T n p ) = nJ(u« T )n p (6) 

where a is the dual of x in Fourier domain and T stands 
for matrix transpose. The above description is schemati- 
cally illustrated in Fig. 3, in which the shaded area repre- 
sents where the power spectrum of an orientated pattern 
is located. 

By substituting Eq. (6) into Eq. (5) we have: 
E{n p ) = | . . . ^J( Ua T )n p \F(u)\ 2 du = n T p Rn p (7) 
where R is a m x m matrix (m — 1 for 2D patterns, and 3 
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Fig. 3. A schematic diagram showing the distribution of the power 
spectrum of an orientated pattern and the derivation of its orientation 
direction. Vector n p shown in the figure is the unit normal vector of a 
Kne that the power spectrum clusters along, whereas d p (n) is the 
shortest Euclidean distance between a point u and the line. 



(8) 



for 3D patterns), given by: 

= j...J(«« T )|F(«)| 2 rf« 

Matrix if will be referred to as the second moment matrix 
of the Fourier spectrum of image f{x). The problem of 
finding the line that the power spectrum clusters along 
is now equivalent to finding n p that minimizes 
function E. 

For simplicity, we shall first look at orientated 2D 
patterns. Denote u-, (i = 1.2) as components of u pro- 
jected onto the two orthogonal- axes; we have: 

R,j = | \u i Uj\F(u [ , u 2 )\ 2 du y du 2 (ij =1,2) (9) 

By applying the power theorem and the fact that differ- 
entiation in the spatial domain corresponds to multi- 
plication by the respective coordinate in the Fourier 
domain [12], Eq. (9) can be simplified and calculated 
directly in the spatial domain [13]; i.e.: 



dx 2 



(10) 



thus, no Fourier transform is involved in the actual cal- 
culation. Furthermore, function E(n p ) defined in Eq. (7) 
is non-negative. Thus: 

njRn p = (n p ,Rn p )^0, V«, (11) 
Eq. (1 1) states that R is positive semi-definite. From Eq. 
(10) it can be shown that R is also real symmetric, thus 
the minimization problem formulated in Eq. (7) is a 
matrix eigenvalue problem and it is solved by n p , the 
vector that is parallel to the eigenvector corresponding 
to the smallest eigenvalue of the second moment matrix 
R. Since R is real symmtric and semi-definite, all eigen- 
values of R are real and non-negative. The smallest 
eigenvalue A min corresponds to the smallest value that 



function E{n p ) can attain, and the eigenvector that cor- 
responds to A min determines the estimated orientation 
direction. 

Similarly, max(£(ji p )) = A majl , where A mM is the maxi- 
mum eigenvalue of R. This suggests the use of the follow- 
ing equation as a measure of anisotropism: 

{t^tf <12) 

For a pattern that is strongly orientated, A raax » A min , 
thus g « 1, whereas for an isotropic pattern A m!U! as Amin 
and g « 0. Analytically, the direction of n p is calculated 

by: 



0(x) =£tan" 



j m 



df\ (df\ 



i((D ! -(l)> 



(13) 



and the anisotropic measure is given by 



{am 


fdf 






<l) (!)"'•}' 


- 


(1IJ 


m 


<£)) 





(14) 

both of which can be calculated directly from the original 
data/(x,y) and its partial derivatives. In Eqs (13) and 
(14) we have substituted x x and x 2 with x and y, respec- 
tively, and n is a local neighbourhood of x = (x,y). 

To demonstrate the effectiveness of the anisotropic 
measure stated by Eq. (14) and how it is used to control 
the shape of the filter kernel, a test pattern shown in 
Fig. 4a has been generated. This pattern has a varying 
anisotropic strength along the x direction. Fig. 4c shows 
the same pattern but with some additive noise. The 
results of applying Eq. (14) to Figs 4a and 4c are given 
in Figs. 4b and 4d, respectively. A fixed window of 7 x 7 
pixels was used for the calculation. It was positioned in 
the vertical centre of the image and shifted horizontally 
from left to right. As can be seen from the figures, the 
anisotropic measure remains high for an image that is 
noise free. Otherwise, g(x) drops as the anisotropic 
strength decreases thus makes the effective support radius 
increase along the horizontal direction. For pixels on the 
left side of the image, a narrower kernel is used to avoid 
destroying image structure during the filtering process. On 
the other hand, a wider kernel is adopted for pixels 
towards the right as the signal varies smoothly, a higher 
efficiency thus can be achieved as the kernel now covers a 
larger area. This adaptiveness to local image structure 
ensures that the efficiency and fidelity conditions men- 
tioned earlier are satisfied throughout the filtering process. 




Z. Yang er at./lmage and Vision Computing 14 (1996) 13S-14S 




Fig. 4. A synthetic orientated pattern (a) with its anisotropic strength (b) measured along the .v axis. The same pattern with additive noise (c) and its 
measured anisotropic strength (d). 



4. Corner recognition 

The anisotropic measure defined above not only gives 
the indication of how strong a pattern is orientated, but 
also provides a convenient way of finding corner and 
junction points within a given image. 

It is known that at edges the anisotropic measure g(x) 
is close to one, thus i>{g(x)) = 1 - g{x) as 0. On the 
other hand, at flat regions the gradient strength |V/(j)| 
is close to zero. It is only at corner and junction points 
that 1 - g{x) » 0 and \ Vf(x)\ » 0. This suggests the use 
of the following expression as a measure of corner strength: 

c(x) = i>(g(x))\mx)\ 2 = (1 - g(x))\Vf(x)\\ (15) 
In Eq. (15), ij>(g{x)) can generally be chosen as a mono- 
tonically decreasing function with respect to g(x) as long 
as it satisfies: 

(16) 

The method stated above involves minimum extra 
computation, and we shall see in the next section that 
our results acquired compare favourably with other 
popular corner detection techniques. For example, one 
well known method relies on the derivative of the 



gradient directions (the curvature of level contours) 
[14] and the equation derived is in the form of: 



+ dx 2 \dy) 



(17) 



To obtain a realistic indication of corner strength, the 
popular Kitchen-Rosenfeld corner detector [15] is 
usually employed which, in effect, is a multiplication of 
the above measure with the local gradient strength, i.e.: 



c{x) 



dy 2 {dxj [dxdy ) Ox dy dx 2 \dy) 
\dx) + \dy) 



(18) 

This involves second order derivatives of the original 
image, and problems may arise if the input data is 
noisy. This is because taking derivatives is a high-pass 
filtering process which amplifies the effect of noise. Other 
available techniques are based on calculating the histo- 
gram of local directional derivatives which entails an 
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analysis of the following different cases: 

• Twin peaks separated by 180° — straight lines or edges; 

• Twin peaks separated by 90° or less — right angle or 
sharp corners; 

• Single peak — termination points for lines; 

• Multiple peaks — junctions. 

Yet again, this method is sensitive to noise and a stable 
result is not always easy to achieve in practice. Further- 
more, it involves a high computational cost as it requires 
histogram peak searching for each local window. 



By incorporating the corner strength c(x) into Eq. (4), 
the proposed adaptive filtering scheme can finally be for- 
malized as: 

<7l(*0) = T 



^(■*o) = (I ~g(x 0 )) — 



(19) 



in which a is a normalization factor that controls how 
faithfully the corners and junctions should be preserved 
during the filtering process. It is usually selected as a 
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Fig. 7. (a) A first order orientated pattern with its Fourier spectrum clustering within a plane that is perpendicular to the orientation direction. The 
shape of the filter is like a thin rod tapering at both ends, (b) A second order oriented pattern with its Fourier spectrum clustering along a line. The 
minimum eigenvalue of the second moment matrix R has a multiplicity of two. The eigenvector corresponding to the maximum eigenvalue is in parallel 
with the line that the power spectrum clusters around. The shape of the filter kernel is this case is like a thin disk which is parallel to the orientation of 
iso-surface layers. 



value between 50% and 200% of the maximum corner 
strength within the image. In Eq. (19) r is the maximum 
support radius, as mentioned earlier. It can be seen that 
as the corner strength c(x) increases, both tr, (jc 0 ) and 
a 2 (x 0 ) decrease to ensure that a small effective window 
is used to avoid smearing. The shape of the kernel, on 
the other hand, is adjusted by the term (1 — g(x 0 )) in the 
expression for ^(jcn). When an edge is encountered the 
term (1 - g(x 0 )) is small. The extent of the kernel is thus 
deformed from a circle to an ellipse, which guarantees 
that the smoothing is performed predominantly along 
but not across the edge. 

5. Results for synthetic images 

The proposed technique has been verified with a syn- 
thetic image. Fig. 5a shows a synthetic image with addi- 
tive noise and features a group of lines with different 
width, a first and a second order grey level strip (the 
increase in pixel intensities along the horizontal direction 
is proportional to x and x 2 , respectively) and five disks 
with different sizes and intensities. All these are 
embedded into a large disk with varying intensity from 
the centre outwards. It is constructed in this way so that 
the results will reflect the filter response with respect to 
different image details. Fig. 5b shows the corner strength 
c(x) when Eq. (15) is applied. As a comparison, Fig. 5c 
gives the result of detected corners by the Kitchen- 
Rosenfeld operator. The contrast between Fig. 5b and 
5c is obvious. Fig. 5c not only has higher background 
noise, but also gives a lot of false responses. Some of the 
lines are also detected as corners because they are 
embedded in a background with varying pixel intensities. 
Our proposed method clearly handles this situation 
much better. As a valuable by-product derived from 



this work, the corner detection scheme mentioned here 
can be a useful tool for other image recognition and 
computer vision applications. 

Figs. 6a-d show the synthetic image with different 
signal-to-noise ratios, i.e. 60 dB, 52 dB, 42 dB and 
28 dB, respectively. The measure of the signal-to-noise 
ratio is defined by: 



s 2 [x)dx 



*?{x)dx 



(20) 



with s(x) and 7(jc) denoting the value of the original 
image and the additive noise, respectively. The input 
image to the filter is f(x) - s(x) + f(x). 

Figs. 6e-h are the processed results of Figs. 6a-d, 
respectively. The thin line at the top of the synthetic 
image has a width of only one pixel, and it can be 
seen that it starts to fade away when the SNR 
decreases. The preservation of edges is generally 
good, though a slight degradation of corners is notice- 
able. The maximum support radius used was 3 pixels, 
and the value a was set to be 75% of the maximum of 
c(x) within the image. An increase of more than 20 dB 
in SNR has been achieved with very minor degrada- 
tion of local features. 



6. Considerations for 3D implementation 

The extension of the technique to 3D requires some 
careful consideration. A 3D image can be considered as a 
2D cine sequence, or a 3D volume data set. There exist 
two situations for 3D orientated patterns: one is that the 
isogrey values of the intensity pattern in 3D space 
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Fig. 8. A single trans-axial head image (a) with its detected corner strength (b) and the filtered image (c). The maximum support radius is 3 pixels, thus 
giving a window size of 7 x 7 pixels. 



constitute parallel lines; and the other is that the isogrey 
values form parallel planes. These two situations are 
illustrated in Figs. 7a and 7b, respectively. The former 
corresponds to conditions in MR imaging where small 
vessels run in parallel within a local area (e.g. peripheral 
blood vessels shown in MR angiograms), whereas the 
latter corresponds to thin surface layers lying in parallel 
with each other. These two cases are termed as the first 
and second order orientated patterns, respectively. For 
first order orientated patterns, the Fourier spectrum 
clusters along a plane through the origin with its normal 
denned by the eigenvector that corresponds to the 
smallest eigenvalue of the second moment matrix R. 
For second order orientated patterns, the Fourier 
spectrum of f(x) shrinks into a line through the origin. 



b 




d 



The smallest eigenvalue of R in this case has a multi- 
plicity of two, thus for any plane that contains the line 
the function E(n p ) will be minimized. The sub-space 
denned by the normal vectors of these planes corre- 
sponds to a plane that is denned by two orthonormal 
eigenvectors of \ min . Since all the eigenvectors of R are 
orthonormal with respect to each other, the eigenvector 
that corresponds to A majl will be parallel to the line 
around which the power spectrum clusters. 

The extension of Eq. (4) to accommodate 3D 
orientated patterns can be subsequently described by 
the following expressions: 

k(x 0 ,x) = / <r- J ^)e-E?..«^)-^^M (21) 
where is the eigenvector that corresponds to X h and: 

ai(Xo) = {' ~ iwfe) }i 1 

(1 = 1,2,3) (22) 

with c(x) as has been denned in Eq. (15). As in Eq. (19), a 
is a normalization factor and r is the maximum support 
radius. They are combined with the corner strength c(x) 
to regulate the effective window size. The three eigen- 
values are used here to control the shape of the filter 
kernel. As shown in Fig. 7, for first order orientated 
patterns, there are one small and two large eigenvalues, 
and the kernel has the shape of a thin rod tapering at 
both ends, whereas for second order orientated patterns 
there are one large and two small eigenvalues; 
consequently, the kernel appears like a disk which is 
parallel to isogrey level surfaces. 



Fig. 9. (a)-(b) Filtered results for the image shown in Fig. 8(a). The 
maximum support radius used for (a) was 4, whereas for (b) was 2. The 
results are very similar to Fig. 8(c). (c)-(d) The differences between 9(a), 
9(b) and Fig. 8(c). The images shown in (c) and (d)have been amplified 
40 times and are displayed with the same brightness and contrast set- 
tings as those for Figs. 9(a) and 9(b). 



7. Application results and discussion 

Several different types of MR images have been used 
to evaluate the performance of the proposed algorithm 
for noise reduction. These images were acquired with a 
0.5T MR scanner (Royal Brompton, London). All 
images in the experiments have the same data array 
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size of 256 x 256 with 1 byte/pixel. They were processed 
on a Silicon Graphics IRIS Indigo R4000 workstation. 

In the first experiment a single trans-axial head image 
from a two-dimensional acquisition is used. The original 
and the filtered images are shown in Figs. 8a and 8c, 
respectively. Fig. 8b is the detected corner strength 
which controls the extent of the filter kernel. The 
improvement of visual quality and structure definition 
is evident. The maximum support radius used in this 
experiment was 3, thus giving a window size of 7 x 7 
pixels. As a comparison, Figs. 9a and 9b provide results 
with r = 2 and r = 4, respectively. The differences 
between Figs. 9a, 9b and Fig. 8c are given in Figs. 9c 
and 9d, respectively, where the subtracted result has been 
amplified 40 times to reveal small discrepancies. It 
appears that the final filtered images do not rely heavily 
on r. This is a desirable feature for images with coarse as 
well as fine structures. The value a for controlling how 
well the corners and junctions are preserved was set to be 
the same as in previous calculations, i.e. 75% of the 
maximum of c{x) within the image. 

Our second experiment used a gradient echo cine 
oblique view of the ascending and descending aorta. Six- 
teen frames were acquired within each cardiac cycle, and 
the 3D version of the algorithm has been applied to the 
data set. Only four images are shown here. They were 
gated at 100 ms, 300 ms, 500 ms and 700 ms after the 
onset of the ECG (electrocardiogram) R wave. The 
filtered images are shown in Figs. 10e~h. These results 
show improved region definition which facilitates the 
region segmentation required for aortic compliance 
studies. 

As the last example, we applied this algorithm to a 
multi-slice spin-echo acquisition of the knee. The images 
are of poor quality due to the thin slice and inferior MR 
receiving coils used. The structures in this data sets have 
many thin muscle surface layers. Using the algorithm 
with the same parameter settings as before, we obtained 
the enhanced images as shown in Fig. 11. There is no 
obvious loss of structural information, and an improve- 
ment of both visual quality and structure definition are 
apparent. 

It is worth commenting on techniques that rely on the 
anisotropic diffusion scheme proposed by Perona and 
Malik [8]. As has been discussed earlier, the method 
allows an image f(x) to evolve over time via the diffusion 
equation: 

= D(x, t)Af{x) 4- VD • (23) 

where D > 0 is a decreasing function with respect to the 
edge strength and A is the Laplacian operator. The result 
is to diffuse f{x) most where the gradient is smallest and 
least where the gradient is largest. Since the orientation 
of the edge has not been taken into account, the efficiency 
of the diffusion process is not ideal. When an input image 



contains a considerable amount of noise, the resulting 
random fluctuations of |V/(jc)| slow down the diffusion 
process. In our experience, we noticed that the filtered 
result can look patchy when very noisy images are 
encountered. Another major drawback of the technique 
is that it distorts sloping edges. As has been shown [16], 
an edge can deform into a multiple or single jagged one 
with a shape that is dependant on selected parameter 
values. Furthermore, the filtered results converge ulti- 
mately to a piecewise constant image definition, thus 
only a small number of iterations should be used if 
subtle image details are to be preserved. The method 
described in this paper, however, does not have all 
these limiting factors. It does not affect the image 
fidelity and can at the same time keep up with the 
efficiency in terms of noise filtering. In addition, all 
parameter values used proved to be robust, and this 
makes the algorithm easily applicable to images with 
varying structures. 



8. Conclusions 

The technique proposed in this paper relies on the 
detection of local orientation features by adaptively 
controlling the filter kernel. The detection of local 
orientation is based on the Fourier analysis of the 
intensity pattern. The evaluation of Fourier trans- 
forms, however, is not necessary for the actual cal- 
culations. A simple relationship between the local 
orientation direction and matrix eigenvectors has 
been obtained. The eigenvalues and their correspond- 
ing eigenvectors are used to control the shape and the 
extent of the filter kernel. Application of this filtering 
process ensures that corners and junctions are well pre- 
served. Experimental results obtained by the applica- 
tion of our method to both synthetic and real images 
show that the algorithm is robust and can cope with 
low signal-to-noise ratio. The noise filtering efficiency 
of the algorithm is good; the computational efficiency, 
however, still needs improvement. The current execu- 
tion time for a 256 x 256 image with 1 byte/pixel is 
about 12-18 seconds, depending on the size of the 
maximum support radius. We are currently implement- 
ing a parallel algorithm to improve the efficiency of our 
method. 
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List of symbols 

x position vector of a point in space 

X; the ith component of x 

u dual of x in Fourier domain 

ttj the rth component of a 

X(, a given point in space 

it a unit normal vector to a line (in 2D) or plane 

(in 3D) 

n p a unit normal vector to a candidate line (in 2D) 
or plane (in 3D) that the power spectrum of an 
orientated pattern clusters along 

u ± a unit vector that is perpendicular to n 

f2 local neighbourhood of x 0 

s( - ) the original image (not corrupted by noise) 

■y( • ) additive noise to s( • ) 

/(•) input image, /(•)=*{■) +7(0 

/,(•) filtered result of /(•) 

F( • ) the Fourier transform of /( • ) 

*(.) filter kernel 

fi( • ) a normalization factor for k( - ) 

p(') a positive and rotationally symmetric cut-off 
function 

r the maximum support radius of k( • ) 

a a normalization factor for adjusting the 

influence of corner strength in k( • ) 
o-j ( • ) a function that controls the major axis of the 

filter kernel 

<r 2 (-) a function that controls the minor axis of the 
filter kernel 

d p (u) the shortest Euclidean distance between a point 
u and a line that has normal vector of n p 

R the second moment matrix of the Fourier 
spectrum of image f{x) 

Rij the (/, ;)th element of matrix R 

A eigenvalue of R 

A, the /'th eigenvalue of R 

A min smallest eigenvalue of R 

A max largest eigenvalue of R 

g{ ■ ) the local anisotropic measure of an orientated 
pattern 

c( • ) corner strength measure 

ip(g(x)) a decreasing function with respect to g(x) used 
for composing c( ■ ) 

■E^iij,) the moment of inertia of the Fourier spectrum 
with respect to a line (or plane in 3D) through 
the origin with a unit normal vector n p 



D diffusion parameter 

T matrix transpose 

A Laplacian operator 

V gradient operator 

| • | magnitude of a complex number or a vector 

d partial derivative operator 
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Abstract 

This paper presents a new operator for corner de- 
tection. This operator uses a variant of the morpholog- 
ical closing operator, which we have called asymmetri- 
cal closing. It consists of the successive application of 
different morphological tranformations using different 
structuring elements. Each of these structuring ele- 
ments used to probe the image under study is tuned to 
affect corners of different orientation and brightness. 
We found that this kind of approach, based on bright- 
ness comparisons, leads to better quality results than 
others and is achieved at a lower computational cost. 

1 Introduction 

Corners constitute attractive 2D features, often 
nsed in computer vision for tasks such as stereovision, 
3D interpretation, motion estimation, and structure 
from motion. They abound in indoor scenes where 
several polyhedral objects and intersecting planes 
(floor, walls, etc.) are present. Corners serve as points 
of interest in two-view matching algorithms [l][2][3]. 
Corner detection is also used in camera calibration for 
the localization of reference points on a calibration 
pattern ([4] for example). 

Corner detection is sometimes realized through the 
analysis of binary edge maps from which chain codes 
are extracted in order to find high curvature points 
[5] [6] [7]. However, most approaches work directly at 
the grayscale level [8]-[17]. These methods usually 
use local measurements in order to obtain a corner 
strength. Non-maxima, suppression and thresholding 
lead then to a binary map showing where corners have 
been detected. These corner finder are usually char- 
acterized by an accuracy of few pixels and a relatively 
high level of false positives. Model-based approaches 
such as [18] [19] also exist and allow corner localization 
at a subpixel accuracy. But these methods are more 
CPU-intensive and are only used after a first corner 
map has been obtained. 

One of the difficulties with corner detection lies in 



the corner definition itself. A restrictive description 
simply defines corners as the junction of two homoge- 
neous regions separated by a high-curvature boundary. 
This definition is incomplete since it does not include 
X, Y and T junctions that should also be categorized 
as corners since they might be the image of 3D corners 
(intersection of three planes). A less rigid definition 
assimilates corners to points with high derivatives in 
several directions. This is a very loose description of 
the term corner since since several "non-corner" points 
fall into this category. 

This paper proposes an approach to corner detec- 
tion based on mathematical morphology. The goal 
was to obtain a fast corner detector that is accurate, 
stable, selective and robust to noise. The next section 
is a short review of existing corner detectors. Section 
3 presents some mathematical morphology concepts. 
Section 4 describes the proposed corner detector and 
Section 5 shows some comparative results. Section 6 
is a conclusion. 

2 Corner detection 

We review here the main corner detectors that work 
directly at the grayscale level. All these methods 
use local measurements in order to obtain a corner 
strength c(x,y) for each point of the image. Local 
non-maxima suppression and thresholding are then 
performed in order to extract points that will be re- 
ported as corners. 

The use of the Hessian determinant of the intensity 
surface to estimate corner strength has been proposed 
by Beaudet [8] . Kitchen and Rosenfeld [9] proposed to 
use the gradient magnitude and the rate of change of 
gradient direction along an edge contour. Very similar 
operators have also been proposed by Dreschler and 
Nagel [10] and Zuniga and Haralick [11]. Deriche and 
Girondon [12] proposed a scale-based approach that 
uses the Beaudet's operator in conjunction with the 
Laplacian. 

Following the idea of points of interest developed 



by Moravec [13], the Plessey detector [14] is based on 
the following matrix: 



M(x,y) = 



(i) 



where (I) denotes a smoothing operation on I. Cor- 
ner strength has been first defined by Noble [15] from 
which a slightly different version has been proposed 
by Harris and Stephen [16]: 

c HS (x,y) = Det(M(x,y)) - k Trace 2 (M(x,y)) 

(2) 

The role of the parameter k is to remove sensitivity to 
strong edges. 

All of the above methods are based on directional 
derivatives. They all suffer from the same drawback: 
local estimation of derivatives is very sensitive to noise 
and, when smoothing is applied, the corner localiza- 
tion precision is reduced. In addition, the computa- 
tional complexity of the smoothing operation, deriva- 
tive estimation and corner strength computation that 
are involved in such methods can be quite high. 

A simpler approach based on brightness compar- 
isons has been proposed by Smith and Brady [17]. 
The SUSAN corner detector is a modified version of 
the edge detector of the same name. It is based on the 
computation of the area of points inside a circular re- 
gion M X y having a brightness similar to the one of the 
center point (a;, y). This area is computed as follows: 



(*.j>eJV„ 



(3) 



The parameter t controls the sensitivity to noise, i.e. it 
defines the similarity between brightness values. The 
value of n(x, y) is therefore compared to a fixed thresh- 
old equal to n max /2 where n max is the maximum value 
that n() can take, that is: 

sv ' J1 [0 otherwise ^ 

The value of this function corresponds to the corner 
strength. In order to reduce the number of false pos- 
itives due to smooth boundary, thin lines and fine 
textures, two criteria must be added. The center of 
gravity of the circular region must be (1) located suf- 
ficiently far away from the center point and (2) all 
pixels lying in a straight line from the center point to 
the center of gravity must be of similar brightness. 



Because of their computational simplicity, methods 
based on brightness comparisons constitute an attrac- 
tive solution to the problem of corner detection. More- 
over, our experiments showed that this kind of ap- 
proach allied accurate corner localization with good 
robustness to noise. The corner detector we propose 
follows this approach and makes use of some basic 
morphological tools. 

3 Mathematical morphology 

Mathematical morphology is a methodology for im- 
age analysis that has been widely used in computer 
vision. The principle of all basic morphological opera- 
tors is to probe the image under study with a structur- 
ing element. This structuring element is a set of pixels 
on which an origin is defined. To evaluate the results 
of a morphological operation on an image point, the 
structuring element is translated in such a way that 
its origin coincides with this image point. The shape 
of the structuring element defines a set Ise^iV) that 
includes all pixels of the image hit by the structuring 
element. Prom this set, the elementary morphologi- 
cal operators erosion and dilation can now be defined. 
The erosion of an image I with a structuring element 
SE is given by: 

Ise(x, y) = min Isefx, y) (5) 
Similarly, the dilation of an image I with a structuring 
element SE is given by: 

Ise& y) = max Ise{x, y) (6) 
Two morphological transformations are defined by 
the successive application of these operators. The 
opening of an image I by a structuring element SE 
is defined as an erosion followed by a dilation: 

where SE is the symmetrical transposition of SE with 
respect to its origin. All image structures that cannot 
contain the structuring element are removed by the 
opening. Therefore, the shape and size of the struc- 
turing element must be set according to the informa- 
tion that is to be extracted. The closing of an image 
I by a structuring element SE is defined as a dilation 
followed by an erosion: 

I c se=(1 5 seYse ( 8 ) 

4 Asymmetrical closing for corner de- 
tection 

In the context of corner detection, one interesting 
choice is to consider a cross-shaped structuring ele- 
ment (Figure 1(a)). Indeed, because of the particu- 
lar shape of this structuring element, the opening and 



closing operators alter mainly this kind of image struc- 
ture. However, these corner detectors suffer from three 
problems: 

1. Opening affects only bright corners over dark 
background while closing affects only dark cor- 
ners over bright background. 

2. Small image structures (including impulsive noise 
and thin lines) are also eliminated and thus can 
be wrongly assimilated to corners. 

3. This kind of corner detection is not rotationally 
invariant. 

A concurrent application of opening and closing can 
solve the first problem, but we introduce a more effi- 
cient solution. We propose to perform what we call 
an asymmetrical closing, that is, the dilation of an 
image using a given structuring element followed by 
an erosion using another structuring element (note 
that asymmetrical opening could also have been con- 
sidered). The central idea is to make dilation and 
erosion complementary in terms of the type of cor- 
ners they affect. This can be realized by choosing a 
cross as the first structuring element and a lozenge for 
the second one (Figure 1(b)). We can then write the 
asymmetrical closing as: 



(9) 



= (4): 



and corner strength will be given by comparing the 
resulting image with the original one, that is: 



c + {I) = |l-I^. i0 | 



(10) 



Basically, the value of c(x, y) corresponds to the 
brightness difference between the corner and its back- 
ground. Note that, as defined, the transformation pro- 
duces a three-pixel L-shaped response in the case of 
dark corners We also observed this kind of multiple- 
pixel response for real images. This must be inter- 
preted as a consequence of the fact that the precise lo- 
cation of smooth corners is not well defined. If needed, 
it is still possible to select only the central point in 
each set. Corners detected by c + on a test image are 
shown in Figure 2(a) (for all experiments to follow, 
we used the structuring elements shown in Figure 1). 
Clearly, rotational invariance and small structure sen- 
sibility have not been solved. To detect the missing 
corners, the following operator can be used (which is 
a 45° rotated version of the preceding one): 



X 



Figure 1: (a) structuring element +. (b) structuring 
element o. (c) structuring element x . (d) structuring 
element □. 



Figure 2(b) shows the corners detected by this oper- 
ator. Surprisingly, it appears that the two comple- 
mentary operators, c + and c x , are sufficient to detect 
corners of almost any orientation. While these two op- 
erators are sensitive to corners of different orientation, 
they are both sensitive to the same small structures. 
Consequently, the combination of these two operators 
should make corner detection almost rotationally in- 
variant and insensitive to small image structures. This 
leads to the following operator: 



c+ , x (i) = |r +i0 -i c x , D | 



(12) 



Cx (i) = |i-r Xi0 | 



(ii) 



Results obtained using this last operator are shown 
in Figure 2(c). This is the operator that we propose 
to use as a corner detector. To extract corners from 
the output of c + , x (I), it appeared to us that a sim- 
ple global thresholding is sufficient even if this leads 
to multiple-pixel response for some corners (i.e. oc- 
currence of a corner represented by a few connected 
pixels in the binary corner map). In fact, we found 
that non-maxima suppression, -which is required in 
the other methods, does not improve the quality of 
the detection in the case of asymmetrical closing. In 
particular, non-maxima suppression does not rule out 
the multiple-pixel response mainly because the cor- 
ner strength at these location are nearly equal (to the 
brightness difference between the corner and its back- 
ground). However, this multiple pixel response be- 
havior is not problematic in most applications and the 
elimination of the non-maxima suppression process re- 
duces the computational load of the corner detection 
task. 

5 Comparative results 

In order to test the validity of the operator c+, x (I) 
as a corner detector, a series of tests were performed. 
Comparisons were made with the Plessey corner de- 
tector and the SUSAN corner detector. 

Sensitivity to noise was tested by adding a gaussian 
noise to the test image. The comparative results are 
shown in Figure 3. These tests demonstrate the su- 
periority of methods based on brightness comparison 
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over differential detectors. This observation has also 
been made in [17]. The SUSAN operator has also been 
reported to be 10 times faster than the Plessey oper- 
ator. We observed that the method proposed here is 
about 2 times more efficient that the SUSAN operator. 

To evaluate the stability of these detectors, we used 
a test sequence of 4 images showing a table from dif- 
ferent points of view. The checkerboard pattern on 
the table creates 24 corners; the goal here is to test 
the ability of each operator to detect these 24 corners. 
Among these corners there are 4 L-junctions, 12 T- 
junctions and 8 X-junctions. All parameters in each 
method have been set in order to obtain the best possi- 
ble results and a comparable density of corner points. 
For each method, the same parameter values are used 
for all images of the sequence. For the SUSAN opera- 
tor and our operator based on asymmetrical closings, 
we have tested two different threshold values. Results 
are shown in Figure 4 where we show 2 of the 4 im- 
ages. Table 1 presents the number of corners that 
each method has detected. It is not surprising to find, 
from the analysis of these results, that a lower thresh- 
old leads to a higher number of detected corners. But, 
at the same time, the number of false positives grows 
rapidly. In order to estimate the number of false pos- 
itives that each method produced, we have counted 
the number of detected corners lying on the table in 
each image. The result of this analysis is presented in 
Table 2. These results demonstrate the reliability of 
the operator based on asymmetrical closings. It ap- 
pears to be more stable while producing fewer false 
positives. It is also interesting to note that, in the 
case of the operator proposed in this paper, several of 
the false positive detections are due to the fact that 
X-junctions tend to produce a two-corner response. 

6 Conclusion 

We have presented a new operator for corner de- 
tection. It uses a variant of the morphological closing 
operator, which we have called asymmetrical closing. 
We found that this kind of approach, based on bright- 
ness comparisons, leads to better quality results than 
other approaches and is achieved at a lower computa- 
tional cost. Stable and accurate corner detection has 
been obtained using the operator presented in this pa- 
per. Because of its algorithmic simplicity, we believe 
that this operator is an efficient means of producing 
input points of interest for feature-based approaches 
to 3D structure and motion estimation problems. 
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Table 1: Number of corners correctly detected on the 
table. A total of 96 corners were considered in the test 
sequence, 24 in each image. 
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(57) ABSTRACT 

A method and arrangement for embedding a watermark in an 
image are disclosed. The watermark consists of a pseudo- 
random, dense subset of image pixels, e.g. a pattern of lines 
(20). Anumber of salient image pixels (21-26), for example, 
local extremes, corners or edges, is identified and it is 
determined whether they lie on (i.e. within a vicinity 5 of) 
the line pattern (21-23) or not (24-26) , In an unwatermarked 
image (FIG. 2A), the number of most salient pixels (21) 
lying on the watermark is substantially the same as the 
number of most salient pixels (25,26) not lying on the 
watermark. The image is watermarked (FIG. 2B) by modi- 
fying the saliency of the salient pixels in such a way that a 
significant majority (21,23) of the most salient pixels (21, 
23,25) is eventually located within the vicinity of the line 



7 Claims, 3 Drawing Sheets 




U.S. Patent 



Feb. 11, 2003 Sheet 1 of 3 



US 6,519,350 Bl 




U.S. Patent 



Feb. 11, 2003 Sheet 2 of 3 



US 6,519,350 Bl 




U.S. Patent Feb. 11, 2003 Sheet 3 of 3 



US 6,519,350 Bl 



3 


4 


3 


2 


9 


3 


4 


3 


5 



2.9 


3.9 


2.9 


1.9 


9.8 


2.9 


3.9 


2.9 


4.9 



-53 



FIG. 5 



3.1 


4.1 


3.1 


2.1 


8.2 


3.1 


4.1 


3.1 


5,1 



Sij=37.8 



H 

■ ■DTj TTjB 
He Kb 



FIG. 6 



us 6,5: 

1 

EMBEDDING WATERMARKS IN IMAGES 

FIELD OF THE INVENTION 

The invention relates to a method of embedding a water- 
mark in an image, comprising the steps of calculating a 
saliency of image pixels, identifying salient image pixels, 
and processing the image in such a way that a predetermined 
percentage of the most salient image pixels lies within the 
vicinity of a predetermined watermark pattern. The inven- 
tion also relates to an arrangement for embedding a water- 
mark in an image. 

BACKGROUND OF THE INVENTION 

A known method of embedding a watermark as defined in 
the opening paragraph is disclosed in M. J. J. B. Maes and 
C. W. A. M. van Overveld: "Digital Watermarking by 
Geometric Warping", Proceedings of the 1998 International 
Conference on Image Processing, Oct. 4-7, 1998, pages 
424-426. In this known method, the watermark is a prede- 
termined image pattern, for example, a pattern of lines. The 
image is watermarked if a statistically high percentage of 
salient pixels of the image lies within the vicinity of the 
watermark pattern. This is achieved by identifying the 
salient pixels, and moving ("warping") them to the vicinity 
of the watermark pattern. The step of image processing thus 
comprises locally changing the geometrical characteristics 
of the image. 

OBJECT AND SUMMARY OF THE INVENTION 
It is an object of the invention to embed the watermark 
using an alternative method of image processing. 

To this end, the method in accordance with the invention 
is characterized in that the step of processing the image 
comprises modifying the saliency of salient pixels. The 
watermark is thus embedded by modifying the saliency of 
image pixels instead of moving them to different positions. 
The saliency is modified by decreasing the saliency of most 
salient pixels not lying within the vicinity of the watermark 
pattern and/or increasing the saliency of salient pixels lying 
within the vicinity of the watermark pattern. Advantageous 
embodiments of identifying and modifying salient pixels are 
defined in the dependent claims. 

BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 shows schematically an arrangement for embed- 
ding a watermark in an image in accordance with the 
invention. 

FIGS. 2A and 2B show watermark patterns and salient 
pixels to illustrate the operation of the arrangement which is 
shown in FIG. 1. 

FIGS. 3Aand3B show histograms to illustrate the opera- 
tion of the arrangement which is shown in FIG. 1. 

FIG. 4 shows a further embodiment of the arrangement 
for embedding a watermark in an image in accordance with 
the invention. 

FIGS. 5 and 6 show sub-images to illustrate the operation 
of alternative options for identifying and modifying salient 

DESCRIPTION OF EMBODIMENTS 
FIG. 1 shows an embodiment of an arrangement for 
embedding a watermark in an image in accordance with the 
invention. The arrangement comprises a salient point extrac- 
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tion module 10, a decision module 11, and an image- 
processing module 12. The arrangement receives an input 
image I and a watermark W, and generates a watermarked 
image l w . 

5 FIG. 2 A shows an example of the watermark. In this 
example, the watermark W is assumed to be a pattern of 
lines 20, but this is not essential. Salient pixels are shown as 
circles 21-26 in the Figure. The diameter of a circle repre- 
sents the saliency of the pixel. A pixel is said to lie in the 

10 vicinity of the watermark, if the distance from that pixel to 
the nearest line 20 is less than a predetermined value. This 
predetermined distance may be a fixed value 8. This defi- 
nition of vicinity, where the lines 20 have a "thickness" 26, 
will be adhered to throughout the rest of the description. The 

I s expression "lie in the vicinity of the watermark" will often 
also be referred to as "lie on the watermark". In FIG. 2A, the 
pixels 21-23 lie on the watermark, whereas pixels 24-26 do 

Alternatively, the expression "lie in the vicinity of the 
20 watermark" can be defined to mean that the saliency-to- 
distance ratio has a predetermined minimum value. Thus, 
pixel 25 in FIG. 2A may lie in the vicinity of the watermark 
whereas, say, pixel 22 does not, because the larger distance 
of pixel 25 to the nearest line 20 is compensated by its larger 
25 saliency. 

The watermark pattern W has a given density. This is 
understood to mean that the watermark pattern including its 
vicinity covers a given portion (p) of the image. Hereinafter, 
it will be assumed that the watermark covers p=50% of the 
image. Accordingly, about p=50% of the salient pixels of an 
unwatermarked image will lie on the watermark. 

The arrangement shown in FIG. 1 operates as follows. 
The salient pixel extraction module 10 comprises saliency- 

35 calculating means 101 for calculating the saliency of the 
image pixels. Embodiments thereof will be described later. 
The module 10 further comprises a selection circuit 102, 
which selects salient pixels. The number of salient pixels is 
small compared with the total number of image pixels. The 

40 salient pixels are reasonably uniformly distributed across the 
image, so as to avoid that modifying the saliency of a pixel 
affects an earlier change of a nearby salient pixel. This can 
be achieved, inter alia, by dividing the image into blocks and 
selecting one or a few salient pixels per block. It can also be 

4J achieved by requiring a minimal distance between salient 
points. As already mentioned, the pixels 21-26 in FIG. 2A 
are salient pixels. 

The salient pixels and their saliencies are applied to the 
decision module 11. This module receives the watermark 

so pattern W and determines which portion of the most salient 
pixels lie on the watermark. In an unwatermarked image, a 
percentage p (here p-50%) of the most salient pixels will he 
on the watermark. In FIG. 2A, the most salient pixels are 
shown as solid circles. Of these most salient pixels, the 

55 pixels 21 lie on the watermark whereas Ihe pixels 25 and 26 
do not. If the percentage of the most salient pixels lying on 
the watermark W does not substantially differ from p=50%, 
the decision module 11 generates a decision signal D=0. In 
response to this signal, the image-processing module 12 

so modifies the saliency of the salient pixels in such a way that 
a significant percentage of the most salient pixels will lie on 
the watermark. Note that the extraction module 10 and 
decision module 11 jointly constitute a watermark detector. 
The image-processing module 12 decreases Ihe saliency 

65 of salient pixels not lying on the watermark and/or increases 
the saliency of salient pixels lying on the watermark. FIG. 
2B shows the saliencies of the salient pixels 21-26 after 
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modification. The saliencies of the pixels 21-23 lying on the 
watermark are increased, the saliencies of the pixels 24-26 
not lying on the watermark are decreased. The process of 
modifying turns "nearly most salient pixels" into most 
salient pixels and vice versa. In FIG. 2B, the pixel 23 lying 
on the watermark is now one of the most salient pixels and 
thus shown as a solid circle. The pixel 26 not lying on the 
watermark is no longer one of the most salient pixels. Thus, 
three of the four most salient pixels lie on the watermark 
after the process of modification. 

FIGS. 3A and 3B illustrate the process in the form of 
histograms. Reference numeral 30 is a graph of the number 
of salient pixels versus saliency. Numeral 31 denotes the 
pixels lying on the watermark pattern W, numeral 32 denotes 
pixels not lying on the watermark. The shaded area denotes 
the set of most salient pixels. In this example, the set 
includes all pixels having a saliency which is larger than a 
given value S r Alternatively, the set may have a predeter- 
mined number of most salient pixels. FIG. 3A shows the 
histogram of an unwatermarked image, having 55% of the 
most salient pixels lying on the watermark. FIG. 3B shows 
the histogram after increasing the saliency of pixels lying on 
the watermark (which shifts line 31 to the right) and decreas- 
ing the saliency of pixels not lying on the watermark (which 
shifts line 32 to the left). After this step, 70% of the most 
salient pixels lies on the watermark. 

In a preferred embodiment of the arrangement, the pro- 
cess of modifying saliency may be repeated until a prede- 
termined majority, e.g. 75%, of the most salient pixels has 
been found to lie on the watermark. Such an embodiment is 
shown in FIG. 4. The arrangement differs from the one 
which is shown in FIG. 1 in that the processed image I w is 
fed back to the extraction module 10 until a significant 
percentage of the most salient pixels has been found to lie on 
the watermark, and the decision module 11 generates the 
signal D=l. 

The step of increasing the saliency of a pixel, which is 
carried out by the image-processing module 12, implies 
locally adding a luminance and/or chrominance image AI to 
the image I in such a way that the saliency S is amplified. 
Similarly, decreasing the saliency of a pixel implies adding 
a luminance and/or chrominance image AI to the image I in 
such a way that the saliency S is attenuated. In view thereof, 
it will be appreciated that the method of modifying is 
strongly related to the method of calculating the saliency. 

In one embodiment of the arrangement in accordance with 
the invention, a 2-dimensional filter forms the saliency- 
calculating means 101. Such a filter will hereinafter be 
represented by a matrix F, for example, the following 3*3 



average of its neighbors. The 3*3 Laplace filter is repre- 
sented by the matrix: 



In this 
1 by adding 



the saliency S, v of a pixel is n 

' J ~. of the matrix F to the image. 



where X is a given weighting factor, I is a 3*3 sub-image 
having the salient pixel in the center, and l m is the modified 
3*3 sub-image. 

FIG. 5 shows an example of this modification process. 
Numeral 51 denotes a 3*3 sub-image with a salient pixel 
having an intensity I lV =9 in the center. The pixel has a 
saliency S—45 in accordance with Eq. 1. Numeral 52 
denotes the sub-image after processing in accordance with 
Eq. 2 and X=0.1, which increases (he saliency to S fJ =52.5. 
Numeral S3 denotes the sub-image after processing in 
accordance with Eq. 3 and X-0.1, which decreases the 
saliency to S- —37.8. 

An alternative method of increasing the saliency S u is 
based on the recognition that S,, is already large, and that the 
sub-image I itself may be used to amplify the saliency, i.e.: 

In another embodiment of the arrangement, the saliency- 
calculating means 101 is formed by an edge and/or corner 
detector, the saliency of a pixel being represented by the 
edge or corner strength. Comer detectors are known per se. 
An advantageous embodiment is described in C. Harris and 
M. Stephens: "A Combined Corner and Edge Detector", 
Proceedings of the 4 rft Aivey Vision Conference, 1988, 
pages 147-151. This corner detector is defined by a matrix: 



F= \ /ft-1 Ao fu 

[/..- A, A-iJ 

The saliency S u of pixel I fJ . (where i and j denote the vertical 55 
and horizontal pixel positions, respectively) is defined by the 
following equation: 



is a Gaussian function with standard deviation a, the symbol 
* denotes convolution, and 

(,m— and 1, = -^ 



S u = /-i,-, + f-yA-i J + /-u l;-i.M + (Eq. 1 ) 

/ft-i hj-i *■ faah.i + /o.i + 

In one embodiment of the invention, the filter F is a Laplace 
filter. This is a high-pass filter which returns values that are 
indicative of local minima and maxima of the pixel values. 
It returns the value zero if the pixel value l u is equal to the 



are the partial derivatives of the image in the directions x and 
60 y, respectively. The matrix M can be relatively easily cal- 
culated by using the following discrete approximations: 
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The matrix M, which can be written in the form 

has a determinant D-AB-C 2 and a trace T-A+B. The comer 
strength R is now defined by: 

where k is a suitable constant, for example, k-0.01. The 
corner strength R is positive for a corner, negative for an 
edge, and approximately zero in a flat region. 

FIG. 6 shows some examples of 5*5 sub-images, the 
center pixels of which were found to have a large corner 
strength. Although the pixels have multi-bit luminance and 
chrominance levels, the sub-images are here shown as 
binary images, that is, pixels having intensities larger than a 
mean or median value are shown in white and pixels having 
intensities less than said mean or median value are shown in 
black. Note that for some sub-images it is immediately clear 
that the center pixel is indeed a comer, whereas for some it 
is not. 

Amplifying the comer strength is achieved by increasing 
the contrast between the pixels that represent the corner and 
the complementary pixels, for example, by adding an 
amount AI to the intensities of the white pixels in FIG, 6 
and/or subtracting an amount AI from the black pixels in 
FIG. 6. Weakening the corner strength is obtained by the 
inverse operation, i.e. subtracting AI from the white pixels 
and adding AI to the black pixels. 

In the foregoing, it has been assumed that an image is 
watermarked if a significant percentage of salient pixels lies 
in the vicinity of the watermark pattern. It will be 
appreciated, however, that the complementary definition 
may be used in practice, i.e. that an image is watermarked 
if a significant percentage of salient pixels lies outside the 
vicinity of the watermark pattern. 

It is further noted that insufficient salient pixels may be 
found in certain areas of the input image. This may be 
particularly the case in uniform areas of synthetic images 
such as cartoons. In this case, salient pixels lying on the 
watermark can be created, inter aha, by adding particular 
noise patterns to said areas of the image, by adding inten- 
sities corresponding to the filter coefficients of the Laplace 
filter, or by artificially creating comers. 

In summary, a method and arrangement for embedding a 
watermark in an image are disclosed. The watermark con- 
sists af a pseudo-random, dense subset of image pixels, e.g. 
a pattern of lines (20). A number of salient image pixels 
(21-26), for example, local extremes, comers or edges, is 
identified and it is determined whether they lie on (i.e . within 
a vicinity 6 of) the line pattern (21-23) or not (24-26). In an 
unwatermarked image (FIG. 2A), the number of mast salient 
pixels (21) lying on the watermark is substantially the same 
as the number of most salient pixels (25,26) not lying on the 
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watermark. The image is watermarked (FIG. 2B) by modi- 
fying the saliency of the salient pixels in such a way that a 
significant majority (21,23) of the most salient pixels (21, 
23,25) is eventually located within the vicinity of the line 
s pattern. 

What is claimed is: 

1. A method of embedding a watermark in an image, 
comprising the steps of calculating (101) a saliency (S iy ) of 
image pixels, identifying (102) salient image pixels (21-26), 

10 and processing the image in such a way that a predetermined 
percentage of the most salient image pixels (21,25,26) lies 
within the vicinity (8) of a predetermined watermark pattern 
(20), wherein said step of image processing includes modi- 
fying (12) the saliency of salient pixels (23,26) by decreas- 

15 ing the saliency of most salient pixels not lying within the 
vicinity of the watermark pattern. 

2. A method of embedding a watermark in an image, 
comprising the steps of calculating (101) a saliency (S ij: ) of 
image pixels, identifying (102) salient image pixels (21-26), 

„ and processing the image in such a way that a predetermined 
percentage of the most salient image pixels (21,25,26) lies 
within the vicinity (5) of a predetermined watermark pattern 
(20), wherein said step of image processing includes modi- 
fying (12) the saliency of salient pixels (23,26), wherein the 
step of calculating the saliency of pixels includes filtering 

25 the image, and the step of image processing includes chang- 
ing the pixel intensities of a sub-image including a salient 
pixel to modify the response of said filter in accordance with 
a desired modification of the saliency of said salient pixel. 

3. A method as claimed in claim 2, wherein the filter is a 
30 2-dimensional filter, and the step of modifying the saliency 

comprises creating a linear combination of the sub-image 
and the filter coefficients of said filter. 

4. A method as claimed in claim 2, wherein the filter is a 
corner detection filter the response of which represents a 

35 corner strength, and the step of modifying comprises chang- 
ing the pixel intensities of the sub-image to modify said 
corner strength. 

5. A method as claimed in claim 1, further comprising the 
step of feeding back the processed image and repeatedly 
carrying out the steps of calculating, identifying and pro- 
cessing until the significant percentage of the most salient 
image pixels lies within the vicinity of the watermark 
pattern. 

6. A method of embedding a watermark in an image, 
comprising the steps of calculating (101) a saliency (S ;J ) of 

4S image pixels, identifying (102) salient image pixels (21-26), 
and processing the image in such a way that a predetermined 
percentage of the most salient image pixels (21,25,26) lies 
within the vicinity (8) of a predetermined watermark pattern 
(20), creating salient pixels in uniform image areas by 

so adding predetermined pixel patterns to said areas, wherein 
said step of image processing includes modifying (12) the 
saliency of salient pixels (23,26). 

7. An arrangement for embedding a watermark (W) in an 
image (I), comprising means (101) for calculating a saliency 

55 (S,j) of image pixels, means (102) for identifying salient 
image pixels, and means (12) for processing the image in 
such a way that a predetermined percentage of the most 
salient image pixels lies within the vicinity (6) of a prede- 
termined watermark pattern (20), wherein said image- 
processing means (12) is arranged to modify the saliency of 

60 salient pixels by decreasing the saliency of most salient 
pixels not lying within the vicinity of the watermark pattern. 



